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ABSTRACT 

We present two-dimensional hydrodynamic simulations of self-gravitating protostellar disks subject 
to axisymmetric, continuing mass loading from an infalling envelope and irradiation from the central 
star, to explore the growth of gravitational instability (GI) and disk fragmentation. We assume that 
the disk is built gradually and smoothly by the infall, resulting in good numerical convergence. We 
confirm that for disks around solar-mass stars, infall at high rates at radii beyond ~ 50 AU leads to 
disk fragmentation. At lower infall rates, however, irradiation suppresses fragmentation. We find that, 
once formed, the fragments or clumps migrate inward on typical type I time scales of ~ 2 x 10 3 yr 
initially, but with a stochastic component superimposed due to their interaction with the Gl-induced 
spiral arms. Migration begins to deviate from the type I time scale when the clump becomes more 
massive than the local disk mass, and/or when the clump begins to form a gap in the disk. As 
they migrate, clumps accrete from the disk at a rate between 10~ 3 to 10 _1 Mjyr _1 , consistent with 
analytic estimates that assume a 1-2 Hill radii cross section. The eventual fates of these clumps, 
however, diverges depending on the migration speed: 3 out of 13 clumps in our simulations become 
massive enough (brown dwarf mass range) to open gaps in the disk and essentially stop migrating; 4 
out of 13 are tidally destroyed during inward migration; 6 out of 13 migrate across the inner boundary 
of the simulated disks. A simple analytic model for clump evolution explains the different fates of 
the clumps and reveals some limitations of 2-D simulations. Overall, our results indicate that fast 
migration, accretion, and tidal destruction of the clumps pose challenges to the scenario of giant planet 
formation by GI in situ, although we cannot address whether or not remnant solid cores can survive 
after tidal stripping. The models where the massive clumps are not disrupted and open gaps may be 
relevant to the formation of close binary systems similar to Kepler 16. A clump formed by Gl-induced 
fragmentation can be as large as 10 AU and as luminous as 2xlO~ 3 L Q , which will be easily detectable 
with ALMA, but its lifetime before dynamically collapsing is only <~ 1000 years. 
Subject headings: accretion disks, planets and satellites: formation, protoplanetary disks, planet-disk 
interactions 



1. INTRODUCTION 

Even a small amount of rotation of a protostellar 
molecular cloud core will result in most of the core mass 
falling onto a disk rather than the central protostar. Un- 
less non-gravitational angular momentum transport is 
very efficient or the initial angular momentum is very 
small, the resulting disk will become quite massive and 
subject to gravitational instability (GI) (Vorobyov & 
Basu 2005, 2006; Rice et al. 2010; Zhu et al. 2010b), 
which can serve as an efficient mechanism to transport 
material inward (e.g., Durisen et al. 2007, and references 
therein). However, GI disks have been shown to be sub- 
ject to fragmentation if the disk cooling time is shorter 
than or comparable to the orbital period (Gammie 2001; 
Rice et al. 2005; Paardekooper et al. 2011). Ignoring ir- 
radiation from the central star/accretion disk, this cool- 
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ing criterion can be satisfied at radii beyond tens of AU 
(Nelson et al. 2000; Boley et al. 2006; Rafikov 2007; Sta- 
matellos & Whitworth 2008; Stamatellos et al. 2010; 
Boley 2009; see Armitage 2010 for a review). 

Gravitational fragmentation has been proposed as a 
mechanism to form giant planets beyond tens of AU 
(Boss 1997; Boley 2009), explaining recent observations 
of the presence of massive planets at large radii (e.g. 
Fomalhaut, Kalas et al. 2008; HR 8799, Marois et al. 
2008) . Alternatively, Gl-induced fragmentation has been 
proposed to explain multiple star systems (Kratter et al. 
2008, 2010a,b; Offner et al. 2010; Hayfield et al. 2011). 
On the other hand, GI fragmentation is complicated by 
the tendency of fragments to migrate inward; if they 
reach the star, they could produce jumps in accretion 
luminosity like those of FU Orionis outbursts (Vorobyov 
& Basu 2005, 2006; but see Zhu et al. 2010a,b for an al- 
ternative model). If the migration stops at the inner disk, 
it may provide a formation mechanism for close binary 
systems (§7.1). 

The variety of clump masses obtained in previous sim- 
ulations suggests Gl-induced fragmentation depends on 
the initial condition of the disks, which should be deter- 
mined by both the rate of infall and the mass of the grow- 
ing central star. A more realistic treatment of GI requires 
consideration of both the way in which mass is loaded 
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onto the disk as infall proceeds, and irradiation heating 
of the outer disk by the central star; the latter is known 
to weaken and/or suppress GI (Cai et al. 2008). Recent 
analytic studies (Levin 2007; Rafikov 2009; Cossins et 
al. 2010; Zhu et al. 2010b) have shown that irradiation 
may even suppress fragmentation beyond 50 AU. How- 
ever, simulations with a large parameter space have not 
been carried out to confirm these analytic studies. 

Investigations of the evolution of clumps have been lim- 
ited so far, largely because of numerical difficulties. Most 
3D hydrodynamic simulations use arithmetically spaced 
computational grids in the radial direction, which can 
only cover a relatively small range of radii. Their limited 
resolution also makes it difficult to resolve the collapsing 
clump. Particle-based SPH simulations can also suffer 
from deficient numbers of particles, especially in 3D (Nel- 
son 2006) . Many initial attempts to treat disk cooling in 
SPH simulations have only applied the orbital cooling 
prescription which leads to a runaway collapse. Recent 
studies have more realistic radiative transfer treatment 
implemented for both grid based (e.g. Boley 2009) and 
SPH (e.g. Stamatellos & Whitworth 2008, Forgan et al. 
2009, 2010, Boley 2009, Rogers & Wadsley 2011) calcu- 
lations. But since 3D simulations are time consuming, 
long timescale high-resolution simulations that can fol- 
low clump evolution are only feasible in 2D. While some 
analytic work has been done (Nayakshin 2010abc) to ex- 
plore long-term clump evolution, the migration process 
is simplified, and further mass growth through gas accre- 
tion is ignored once the clump has formed. Furthermore, 
2D simulations enable us to study clump formation in a 
large parameter space of both disk and infall properties. 

In order to study Gl-induced disk fragmentation under 
the influence of both irradiation and mass infall, we have 
carried out high-resolution two-dimensional (2D) hydro- 
dynamic simulations with thermal physics included sim- 
ply. We show how fragmentation depends upon the mass 
infall rate and the mass infall radius. Furthermore, with 
high resolution long-timescale simulations, we are able to 
trace the clump evolution. We find that clumps migrate 
inward rapidly, unless they become massive enough to 
open gaps, and the accretion rate onto the clump is usu- 
ally high due to the large cross section. Inward migrating 
clumps can be subject to tidal destruction as they mi- 
grate inward because of their shrinking Hill radii. But 
the detailed accretion physics around the clump and the 
dynamically unstable second core phase complicates this 
issue. 

In §2, we introduce our numerical model. In §3 and 
4, we present our simulation results regarding disk frag- 
mentation and clump evolution, in §5 we discuss the fate 
of the clumps, while in §6, we develop analytic models 
for disk fragmentation and clump evolution to compare 
with numerical results. The implications are discussed 
in §7, and conclusions are drawn in §8. 

2. MODELS 

We use the FARGO- ADSG code (Baruteau & Mas- 
set 2008) in a two dimensional (R, <j>) fixed polar grid. 
FARGO- ADSG is an extended version of FARGO (Mas- 
set 2000) , which uses a finite difference scheme with stan- 
dard source and transport steps that are similar to ZEUS 
(Stone et al. 1992). However, additional orbital advec- 
tion has been applied which removes the average az- 



imuthal velocity so that truncation error is reduced and 
the timestep allowed by the Courant-Fricdrichs-Lewy 
condition is significantly increased. Thus FARGO- ADSG 
enables us to perform high resolution simulations over 
long timescales. Furthermore, the effect of the central 
star acceleration due to the non-axisymmetric potential 
of the disk is included as an indirect potential term in 
the fluid equations (Baruteau & Masset 2008). 

2.1. Radiative cooling and infall 

Starting from the publicly available version of FARGO- 
ADSG, we have implemented the hydrodynamic equa- 
tions with both a simplified radiative cooling treatment 
using detailed opacities (Appendix A, updated from Zhu 
et al. 2009b) and mass infall. The hydrodynamic equa- 
tions are 
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where £ is the disk surface density, E = J edz is the 
internal energy per unit area, P = J pdz is the verti- 
cally integrated pressure, n is the viscous stress tensor 
(which is turned off in this paper) , and v and <E> are the 
velocity and gravitational potential (including the self- 
gravity potential), v is constant along z and $ is de- 
rived from the thin disk approximation. An ideal gas 
equation of state with 7=7/5 has been assumed, so that 
P = E(*y — 1). Then the averaged disk 2D temperature 
is T 2 d = Pm/=KE, where 3? is the gas constant. It can 
be shown that if T{Z) A cx Tq J pdZ in the vertical di- 
rection of a 3D disk, which is close to the structure of a 
viscous heating dominated disk with a constant opacity, 
the disk's central temperature T c = 4/5T 2 _d where T 2 d 
is calculated above with the total pressure and surface 
density. With irradiation included, T 2 jj is closer to T c , 
so we set T 2 d = T c . f; n , Fj n and E; n model the effects of 
infall on the disk mass, angular momentum, and energy. 
These terms are described below in Eqs. 6-8. 
The cooling rate per unit area, Ac, is modeled via 
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where r = (Y,/2)kr(p c ,T c ) is the optical depth to the 
disk midplane at radius R. The Rosseland mean opac- 
ity k(p c ,T c ) is calculated using the opacity table given 
in Appendix A with the disk midplane density p c =S/2H 
and the midplane temperature T c . The particular form 
t/(1 + t 2 ) is chosen so that the cooling term has the cor- 
rect form in the optically-thick and -thin limits. (In the 
optically thin limit, n should be the Planck mean opac- 
ity. Because the dominant opacity is from dust, which 
has a relatively slow variation with wavelength, we just 
use the Rosseland mean opacity for simplicity.) 

The term aT^ Kt represents the heating effect of the ir- 
radiation from the central star, assumed to vary as 
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where L* is the total luminosity of the star and f(R) ac- 
counts for the normal component of the irradiation from 
the central star to the disk. Here we set f(R) to have 
the constant value f(R) — 0.1, which is roughly what 
is expected for flared disks (Kcnyon & Hartmann 1987; 
Chiang & Goldreich 1997; D'Alessio et al. 1998), and 
=L©. For disks without internal heating by spiral 
shocks, the temperature will relax to T ext based on Eq. 
4. 

The effect of infall onto the disk is represented by fj n , 
F in and E in in Eqs. 1-3. f; n and E in are the rate of disk 
mass and internal energy increase per unit area due to 
the infall, which can be directly derived from an envelope 
infall model. The form of F; n , however, is less immedi- 
ately apparent; it is the equivalent shear force induced 
from merging infalling mass with the disk mass when 
they have different velocities. The relationship between 
F in and the azimuthal velocity of the infalling envelope 
can be derived if we write the angular momentum equa- 
tion from Eqs. 1 and 2 
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where J = T,Rvg. Thus R(£ m vg + F- lnt g) is the additional 
angular momentum brought by the envelope falling onto 
the disk per unit area, which should be Rfi n Vi n ,6, where 
v\ nt g is the azimuthal velocity of the infalling envelope as 
it lands on the disk at radius R. If the infalling enve- 
lope has the same velocity as the disk rotational velocity 
{vin,e = vg), there is no shear between the disk and the 
infalling envelope and Fi nj e=0. 

In our simulations, we minimize the envelope's impact 
on the disk by assuming F^g = in order to study the 
disk's response just to the mass loading from the infalling 
envelope. This is actually not a bad approximation based 
on the rigid rotating envelope infall model of Terebey et 
al. (1984) and Cassen & Moosman (1981), which suggest 
most of the infalling material falls onto the centrifugal 
radius (R c ) where it has just slightly less than Keplerian 
rotation (i>i n ,0 ~ vg). This also implies F; n .fl=0. In our 
simulation we add mass to the disk within a small range 
in radii (R a to Rb) according to 



fin(-R) 



M in (1- P )R-p 
2ttR rI~p - r\-p 



(7) 



so that by integrating over R a to Rb the total infall rate 
is Mj n . The parameter p is chosen to be 0.75 so that the 
disk within R a and Rb has the same Q during the infall 
as if the disk is irradiation dominated (see Eq. 5). We 
also assume that the infalling material has a temperature 
equal to T cxt such that 

E in (fl) = f in (fl) ^ . (8) 

(7 - 1)/ZTO H 

This is not precisely correct; backwarming from the enve- 
lope can actually heat the disk. We ignore this effect for 
simplicity in our simulations, but we discuss its possible 
impact in §7.2. 

2.2. Initial conditions and simulation setup 



Currently, most simulations studying GI start with 
massive and non-irradiated disks. In such cases the oc- 
curence or otherwise of fragmentation depends upon the 
initial conditions (Boley 2009). Without external heat- 
ing (irradiation), a disk with a finite initial tempera- 
ture can cool to a point such that the Toomre param- 
eter Q = c s k/ttGT, ~ 1, where k is the epicyclic fre- 
quency, and the disk becomes gravitationally unstable. 
The equivalent orbital cooling time (t coo \) sensitively de- 
pends on the initial surface density distribution, 
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In the optically thick, viscous heating dominated case 
= (3/8)T^ ff r and the cooling time becomes 

nsc„r c 3scr e n 1 

<W" ~ CT 16T C 4 a S4-2aJj.10.5-3a ' ^ W > 

where we have assumed that Q is constant to constrain 
the temperature , kr — CT®, where C is a constant, 
and the disk follows the Keplerian rotation law (so that 
k = O = CIk) in obtaining the final proportionality re- 
lation. Since a = 1.5 for our grain opacity (Appendix 
A), larger initial surface densities at larger radii lead to 
shorter cooling times and fragmentation. 

During a realistic star formation process the disk is 
built up over a finite time from the infall of the parent 
molecular cloud core. Thus the initial conditions depend 
on the mass and angular momentum distribution within 
the cloud, which can be quite complicated, and in any 
event are poorly constrained. Moreover, the thermal his- 
tory of the disk can affect disk fragmentation (Clarke et 
al. 2007). To simulate the disk mass growth from infall 
until it becomes gravitationally unstable we start with 
a relatively low mass disk and then gradually add mass 
to the disk at a constant rate. When the disk becomes 
massive enough and gravitationally unstable, spiral arms 
generated by the GI will transport the infalling mass to 
the central star at a rate matching the infall rate if the 
system can achieve a steady-state. With this procedure, 
we avoid strong relaxation at the beginning of our simu- 
lations which may cause fragmentation due only to initial 
conditions; this turns out to be important for simulation 
convergence, discussed in §2.3. It also allows us to control 
the infall rates and infall radii as two free parameters. 

We have performed a variety of simulations with infall 
at different rates and radii. Since the disk undergoing 
GI tries to accrete to match the infall rate, M <~ M- m 
in order to maintain Q ~1, we can vary M; n to control 
the disk accretion rate. These simulations are summa- 
rized in Table 1. Different cases are named as R+infall 
radii+infall rates. Thus R40_3e-4 means the typical in- 
fall radii is 40 AU, while the infall rate at <~ 40 AU is 
3 x 10~ 4 M Q yr" 1 . The special cases have notations after 
the infall rates which are self-explanatory or explained in 
the footnote. In the standard cases, the azimuthal grid 
resolution is 512. We use logarithmically spaced radial 
grids and choose the radial resolution such that every 
grid cell has the same azimuthal and radial size. In most 
of our simulations the radial resolution is ~ 400, with the 
inner radius at 5 or 10 AU and the outer radius at 1000 or 
1500 AU (see Table 1). We set /(i?)=0.1 in Eq. 5, so that 
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H/R ~ 0.029( J R/AU) 25 (or T irr (l AU)=215 K with ir- 
radiation heating only). The initial disk surface density 
is E = 4.16 x 10 4 (i?/AU)' 175 g cm~ 2 with an exponential 
cut-off beyond Rb (obtained by multiplying the surface 
density by the factor exp(5 — 5R/Rb) if R > Rb, where 
Rb is defined in Eq. 7). With this setup, the disk is grav- 
itationally stable (Q ^2) within Rb, and the disk surface 
density is extremely small close to the outer boundary 
at <~ 1000 AU. With infall, the disk gradually becomes 
gravitationally unstable near infall radii; spiral arms de- 
velop and start transporting mass inwards, resulting in 
the inner disk gradually becoming gravitationally unsta- 
ble. 

We ran our simulations for up to 5xl0 4 years for non- 
fragmenting cases. If the disk fragments, we run the sim- 
ulation until a clump is close to the inner boundary and 
causes negative densities (details in Appendix B). (In the 
R200_3e-6 case we ran the simulation for 2xl0 5 years to 
test whether this disk fragments on a longer timescale.) 
A timescale of 5 x 10 4 years is not much smaller than the 
estimated length of the infall phase for low-mass pro- 
tostars, and it is equivalent to 17 orbits at our largest 
infall radius (200 AU), long enough to see if the disk will 
fragment. 

We assume the central protostar has a fixed mass (M Q ) 
throughout our simulations in order to study the disk re- 
sponse to the mass infall only. Thus our simulations are 
approximating the disk evolution near the end of proto- 
stellar collapse, when the central object is quite massive 
and emits significant radiation. The expected mass of 
the GI disk from i?; n to R out can be estimated by assum- 
ing Q — 1 with the temperature set by our irradiation 
prescription (Eq. 5): 

M - a22 ((w)° 25 -(^)°' 2S ) M - 

With R out =100 AU and R in =10 AU, we have M d = 0.3 
M . Thus the disk is massive compared with the central 
star if irradiation is considered. 

To study how irradiation affects fragmentation, we ran 
three additional simulations with high, low, and no irra- 
diation, respectively. The high irradiation case, R200_3e- 
5L100, is similar to R200_3e-5 but with 100 Lq, typical 
of an A type star. On the contrary, the decreasing irradi- 
ation case, R100_le-5I, is similar to case R100_le-5, but 
after 5xl0 4 years, with the same infall rate, the irradia- 
tion parameter in Eq. 5 gradually decreases at a rate 

f(R)= 0.1 (l \ ^ , (12) 

J V ' \ 5 x 10 4 years J ' V ' 

until the disk fragments. Equation 12 is chosen so that 
the disk scale height due to irradiation decreases linearly 
with time and the disk mass decreases linearly with time 
(assuming Q =constant). The resulting mass accretion 
rate due to the decreasing irradiation is much smaller 
than the infall rate at R>85 AU, so the disk accretion 
rate is still controlled by the infall rate. Finally, we per- 
formed one simulation similar to R200_3e-6 but without 
irradiation at all, which is labeled as R200_3e-6noirr. The 
entire suite of parameters explored in our simulations is 
summarized in Table 1. 



2.3. Numerical tests 

To explore how robust our simulations are, we have 
carried out simulations with different settings. The most 
important issue is the convergence on the fragmentation 
criterion raised by Meru & Bate 2011. Some other issues 
include the 7 value and the gravitational force smoothing 
length. 

2.3.1. Resolution tests 

To study whether the critical fragmentation radius de- 
pends on the numerical resolution, we have re-calculated 
several marginal fragmentation and non-fragmentation 
cases, but with double the resolution in both radial 
and azimuthal directions as shown in Table 1. Non- 
fragmentation cases remain non-fragmenting at higher 
resolution and the fragmentation case still fragments. A 
particular concern applies to the non-fragmentation case, 
in which numerical dissipation may prevent disk frag- 
mentation. Thus, for R100_le-5 case, we double and even 
quadruple the resolution in both radial and azimuthal di- 
rection, but the disk still remains non-fragmenting. We 
want to point out that the irradiated disk has larger scale 
height which is better resolved numerically. We conclude 
that an azimuthal resolution of 512 cells is enough to 
study disk fragmentation, which is expected considering 
each disk scale height is resolved by 7.5 grid cells at 100 
AU (considering H/R ~ 0.029(i?/AU) - 25 ) and the neu- 
tral stability critical wavelength (L = 2nH for Q = 1) is 
resolved by 45 grid cells, far larger than that required by 
the condition given in Nelson (2006). 

Recent SPH simulations by Meru & Bate (2011) have 
questioned the convergence of the orbital cooling frag- 
mentation criterion for GI disks. However, in our resolu- 
tion tests from 512 to 2048 azimuthal grid cells and the 
corresponding quadruple resolution increase in the radial 
direction, our results converged for both non-fragmenting 
and fragmenting disks. We attribute this discrepancy to 
the fact that we gradually build the disk from a grav- 
itationally stable state to an unstable state by adding 
mass to the disk. This process not only simulates the 
real astronomical condition of infall but also allows the 
disk to regulate itself. In a sense this is similar to Clark 
et al. (2007), who increased disk cooling gradually. Sim- 
ulations by Paardekooper et al. (2011) confirmed that 
the treatment in Clark et al. (2007) leads to a conver- 
gent cooling criterion, which is consistent with our find- 
ing that gradually built-up disks converge. 

2.3.2. Adiabatic index 

The two dimensional adiabatic index 7 is related to the 
three dimensional T in the static limit (e.g., Goldreich et 
al. 1986). In a strongly self-gravitating disk, 7 = 3 — 2/T. 
If the three dimensional V = 1.4, the two dimensional 
7=ll/7<~ 1.57. We explore both the quasi-steady GI 
and highly dynamic fragmenting states, so it is not obvi- 
ous what 7 value we should choose. We have simulated 
two cases, non-fragmenting and fragmenting disks, with 
both 7=1.4 and 1.57 (Table 1). Our results suggest that 
the fragmentation criterion does not depend sensitively 
on 7 in the range from 1.4 to 1.57. This is consistent 
with Rice et al. (2005). On the other hand, Boley et al. 
(2007) have found that variation of 7 due to thermaliza- 
tion of the rotational levels of molecular hydrogen can 
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affect gravitational instability (Boley et al. 2007). In our 
already simplified 2D simulations, however, we simply 
assume 7 is constant. 

2.3.3. Gravitational force smoothing length 

To avoid small scale density fluctuations leading to a 
run-away collapse, we use a smoothed gravitational force 
(Baruteau & Masset 2008). Considering that the disk 
vertical structure can smooth the disk's global gravita- 
tional potential, we assume a smoothing length of 0.0264 
R, which is 0.6 disk scale height at 5 AU and 0.3 disk 
scale height at 100 AU. Since a large smoothing length 
may lead to bigger clumps, we have carried out one sim- 
ulation with smoothing length of 10 times smaller. The 
clump mass from this simulation turns out to be close to 
the clump mass from the corresponding large smoothing 
length simulation. This is consistent with Baruteau et 
al. 's (2011) finding that a small smoothing length does 
not change the properties of the gravo-turbulence. 

3. RESULTS: FRAGMENTATION 

Images of the disk surface density for models with 
different infall rates and infall radii at intermediate 
times are shown in Figure 1. In this figure, the infall 
rate increases horizontally from 3xlO _6 M yr _1 (left) to 
3x 1O _4 M0 yr _1 (right). The infall radius runs vertically 
from 65 AU (top) to 200 AU (bottom). These runs are 
listed from R65_3e-4 to R200_e3e-6 in Table 1. At a given 
infall rate (moving from top to bottom down a column in 
Fig. 1), the spiral arms become narrower and more unsta- 
ble with increasing infall radius and eventually break into 
fragments. At marginal fragmentation only one clump 
forms, while for a more unstable case, several clumps 
form simultaneously. Higher infall rates at the same ra- 
dius (moving from left to right along a row in Fig. 1) also 
lead to more unstable spiral arms. Similarly, fragmenta- 
tion becomes more severe (one clump to multiple clumps) 
with higher infall rates at the same radius. The critical 
radius where the disk starts to fragment increases with 
decreasing infall rate (marginal fragmentation cases form 
a diagonal line in Fig. 1 moving from top right toward 
bottom left). 

Figure 2 summarizes our fragmentation results, and 
they are in approximate agreement with our previous an- 
alytic estimates for the critical radius for fragmentation 
as a function of the mass infall rate (Zhu et al. 2010b, 
and further details in §6). These estimates are shown as 
curves in Fig. 2. 

The special case R100_le-5I (decreasing irradiation) 
is shown in Fig. 3. This disk fragments at R=400 AU 
when f(R) decreases to 0.013. A comparison simulation 
that maintains full irradiation (run R100_le-5), which is 
shown in the lower right panel of Figure 3, does not frag- 
ment. We thus confirm that irradiation tends to suppress 
fragmentation at fixed infall rate. In the zero irradiation 
limit model R200_3e-6noirr fragments, unlike its coun- 
terpart with full irradiation, which again confirms that 
irradiation suppresses the fragmentation. 

In the high irradiation limit (run R200_3e-5L100) the 
disk did not fragment by the end of our simulation 
(5xl0 4 years), unlike its 1 L Q counterpart. Again: irra- 
diation suppresses fragmentation. This also suggests that 
for high- luminosity systems (where L ~ 100 L Q ) gravi- 
tational instability can transfer mass inwards smoothly 



from hundreds of AU; an analytic treatment of this result 
is given in §6.1. 

4. RESULTS: CLUMP EVOLUTION 

In a fragmenting disk, clump formation and evolution 
should resemble protostar formation and evolution (Lar- 
son 1969, Bodenheimer et al. 1980). During disk/cloud 
fragmentation, a quasi-hydrostatic core forms (this stage 
is called the "first core" in protostellar studies). When 
the core's central temperature rises from ~1500-2000 K, 
the dust sublimates, then hydrogen molecules begin to 
dissociate. At this point the ratio of specific heats 7 de- 
creases below 4/3, the clump becomes dynamically un- 
stable and it collapses to form the so-called "second core" 
(Larson 1969). 

In our 2D simulations we cannot simulate the dynam- 
ically unstable "second core", and whenever dust sub- 
limates at the midplane our assumption that the disk 
has uniform opacity vertically becomes invalid. Thus, 
although we allow the clump to be heated up above 1500 
K, our results at those stages are not accurate. Whenever 
the clump is above 2000 K and the center of the clump 
has collapsed to the "second core" , a circum-clump disk 
will form, and its accretion becomes important for the 
clump's further growth. In this paper, we focus on the 
"first core" stage. Even for this stage, we found there are 
three distinct fates for the clumps. During the clump's 
inward migration, if the clump grows massive enough to 
open a gap, it slows down or even stops migration. On 
the other hand, if the clump fails to open a gap, contin- 
uous migration can eventually lead to tidal destruction. 
There is a third possibility that the clump manages to 
migrate inward all the way to the central star without be- 
ing tidally destroyed. Since our inner boundary is at <~ 10 
AU, we can not follow the clump all the way to the cen- 
tral star to confirm this possibility, but provide further 
analytic estimates to investigate this in §6.2. Nonethe- 
less, our simulations do exhibit clumps that manage to 
reach the inner boundary of the computational domain 
without being tidally disrupted. 

We want to emphasize that all the clumps studied here- 
inafter are limited to the marginally fragmentation cases, 
in which only one or two clumps form in the disk. They 
do not suffer the strong clump-clump scattering seen in 
cases of multiple fragmentation, which complicates their 
evolution. 

4.1. Clump structure 

Figure 4 shows the clump from the run R100_3e-5 with 
the numerical grid overlaid on the surface density dis- 
tribution. The left panel shows the clump when it first 
forms at 150 AU after 2.2x 10 4 years, and the right panel 
shows its density 2000 years later when it has migrated 
in to 50 AU. The clump is well resolved at both times. 

Figure 5 shows the temperature distribution of the 
clump (corresponding to the right panel of Figure 4). 
Velocity vectors are superimposed showing the velocity 
field measured in a frame rotating at the same rate as 
the centre of the clump. Shocked regions form at the 
edges of the clump where the infalling material from the 
disk collides with the pressure-supported clump (a sim- 
ilar phenomenon is observed in the ID simulations by 
Larson 1969). As the clump accretes gas with high rela- 
tive angular momentum, it develops significant rotation 
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(2 km s -1 , compared with the Keplerian rotation speed 
of 4 km s _1 in the clump rotating reference frame, as 
discussed in detail below). The prograde clump rota- 
tion is similar to that seen in simulations of circumplan- 
etary disks (Lubow ct al. 1999), presumably because of 
a similar flow pattern occurring as gas enters the clump 
Hill sphere from the surrounding protostellar disk. This 
is also consistent with the results of Voroybov & Basu 
(2010) and Boley et al. (2010). 

The clump's structure at the different times shown in 
Fig. 4 are revealed in detail in Fig. 6. The clump's cen- 
tral density and temperature increase with time during 
the collapse of the clump. For a hydrostatic core, the 
clump's self-gravity should be balanced by the sum of 
the pressure and centrifugal forces due to the clump's 
rotation. When the clump first collapses, it collapses ra- 
dially until it is pressure supported. Since at this early 
stage little mass has been accreted that adds angular 
momentum to the clump, the clump rotates slowly. The 
pressure gradient dP/dR then balances the clump's self- 
gravity and is almost one order of magnitude larger than 
the centrifugal force (left panel of Fig. 4). Later, loss of 
thermal energy due to radiative cooling leads to further 
clump contraction, and the clump spins up. Disk mass 
with high specific angular momentum is also accreted, in- 
creasing the rotation of the clump, and centrifugal force 
becomes comparable to the pressure gradient as shown 
in the right panels of Fig. 6. In order to qualitatively 
estimate how much angular momentum is brought by 
the accreted gas from the disk, an analogy can be drawn 
between accreting giant planets in disk gaps and our ac- 
creting clumps. Ayliffe & Bate (2009) have shown that 
the accreted gas from the disk falling onto the planet 
forms a circumplanetary disk with a radius 1/3 of the 
planet's Hill radius. Thus if the angular momentum of 
the accreted gas dominates our clump's initial angular 
momentum, our clump is expected to be close to or even 
larger than 1/3 of the planet's Hill radius considering 
the clump is also partly pressure supported. This is con- 
firmed in our simulations, shown below, where our clump 
radii are close to half their Hill radii. 

To study clump evolution further, we define a clump 
mass, M c , and radius, r c , as follows. First, wc choose 
a large circle around the clump and integrate the total 
mass. If the Hill radius r H = i?(M c /3M st ) 1 / 3 calculated 
by this mass is larger than the the initial radius of the 
circle, we choose a smaller circle by iteration until the 
Hill radius is the same as the radius of the circle. Then 
we define the clump center as the position of peak surface 
density and we compute the average disk surface density 
around the circumference of the clump's Hill radius. 
After this step, we look for the closed curve r(9) around 
the clump center where its surface density is equal to 
0.5 * (\og w Y, H + \og 10 ^ max ), where T, max is the sur- 
face density at the clump center. Wc then define the 
clump radius as the average distance from the curve to 
the clump center (r c —J r{9)d6 /2tt , where 9 is the angle 
between the line connecting each position on the curve 
and the clump center with respect to a reference direc- 
tion, taken to be the line joining the the central star to 
the clump) . This clump radius can be used to determine 
whether or not the clump can be tidally destroyed, by 
comparing it with the Hill radius, which will be discussed 



in detail in §5.2. 

Quantitatively, the clump in the left panel of Fig. 4 for 
R100_3e-5 case is at 139 AU with mass 0.026 M and Hill 
radius 26 AU. The total specific angular momentum J 
with respect to the clump center is 4x 10 17 cm 2 s" 1 which 
is similar to Boley et al. (2010). When the clump accretes 
mass, both mass and angular momentum increase. For 
the clump in the right panel of Fig. 4, its mass is 0.159 
M (the upper panel in Fig. 7) and the specific angular 
momentum is 10 19 cm 2 s _1 . Such high specific angular 
momentum can limit contraction of the clump even in the 
"first core" stage. For an order of magnitude estimate, 
the clump in the right panel cannot contract to become 
smaller than ^AU scale. This is roughly consistent with 
the 8 AU clump radius calculated using the procedure 
defined in the last paragraph. Even with the significant 
growth in mass, the core radius remains approximately 
half the Hill radius at all times, as shown by the right- 
hand vertical axis in the middle panel of Fig. 7. However, 
the angular momentum distribution within a clump and 
how it is transported outwards is unclear. We do not 
include explcit viscosity in these simulations, so evolution 
of the clump spin angular momentum distribution can 
arise through tides and numerical diffusion. 

Considering that the clump's dynamical timescale 
(free-fall timescale) is significantly shorter than its 
Kelvin-Helmholtz timescale (Nayakshin 2010a), the 
clump can be approximated as quasi-hydrostatic during 
its collapse (more timescale estimates are discussed in 
§5.2). This is confirmed by the bottom panel of Fig. 7 
where GM c /^Rr c T c ~ 1.5 is almost constant throughout 
the clump evolution. As shown in Eq. C7 in Appendix 
C, for a spherical, hydrostatic, radiative clump, we ex- 
pect GM c /$tr c T c <~ 0.7. However, considering the rather 
rough definition of r c in our 2D simulation, the geometric 
difference between 2D and 3D, and using the disk mid- 
plane opacity to represent the envelope opacity at that 
disk radius, a factor of two uncertainty is expected and 
can be roughly considered as our error bar. 

4.2. Clump accretion 

In the previous subsection we explored clump struc- 
ture by following an individual clump. From this sub- 
section forward we will study general clump evolution 
(e.g. accretion and migration) . After the disk starts frag- 
menting, a clump quickly evolves to a quasi-hydrostatic 
state. As shown in the previous section, the core radius 
is approximately half of the Hill radius, leaving plenty 
of space for disk material to accrete onto the collapsed 
clump. Due to the stochastic nature of the clump accre- 
tion and migration in GI disks, we will follow the clumps 
from a number of disk fragmentation cases to increase 
our sample. 

Clump accretion rates as a function of time are shown 
in the left panel of Fig. 8. They vary from 10 _5 M Q yr^ 1 
at 200 AU to 10- 4 M Q yr- 1 at tens of AU. These accre- 
tion rates can be estimated by assuming that the accre- 
tion cross section corresponds to 1-2 clump Hill radii: 

M c = 2E / dxv shcar (13) 
Jo 

where u sho ar = 3/2fix, x = R - i? c ium P and i? c i ump is the 
clump's radial position in the disk. a; max corresponds to 
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the impact parameter where the last unbound orbit sits. 
This impact parameter has been studied in the context 
of planetary rings by Petit & Henon (1986) and Henon & 
Petit (1986), and in giant molecular cloud scattering by 
Gammie et al. (1991). We can understand the gas flow in 
a gaseous Keplerian rotating disk by applying the results 
from the circular restricted three-body problem. If disk 
material is orbiting at a radius very close to the clump 
(small impact parameter), the Jacobi integral is large and 
the disk material is restricted to the horseshoe region in 
the vicinity of the corotation radius. If disk material 
orbits far from the clump, it will follow its Keplerian 
orbit with very little perturbation from the gravity of 
the clump. Only disk material with intermediate impact 
parameters can be captured. The impact parameter of 
the last unbound orbit as in Eq. 13 normally lies between 
1-2 Hill radii (Gammie et al. 1991; Bryden et al. 1999; 
Lubow et al. 1999). The planet-disk interaction studies of 
Bryden et al. (1999) and Lubow et al. (1999) have clearly 
demonstrated how the planet accretes from the gaseous 
disk. Although spiral shocks complicate the orbits of the 
disk material, gas in the circumstellar disk orbiting 1-2 
Hill radii away from the planet is deflected into the Hill 
sphere by the planet's gravity. In a disk undergoing GI, 
the clump forms within a high density spiral arm, and 
mass can also flow along this arm onto the clump, as 
indicated by Fig. 5. Even with this further complication, 
we find that Eq. 13 can describe the accretion rate onto 
a clump accurately. 
Assuming x max = f c rn, we can estimate 



the expression (Tanaka et al. 2002) 

R h 2 2tt 



R ^Cqn ft 



(16) 



where C = 3.2 + 1.468£, £ is the slope of the surface 
density S oc r" ? , h = H/R, fj, = ttY<(R)R 2 /M„, and 
q = M c /M*, where M c is the clump mass. With our disk 
parameters, this becomes 



T~mia — 784 



0.01M, 







R 



100 AU 



1.75 



yrs. 



(17) 



M c = -Y.n{f c r H f 



(14) 



where higher order terms in x/R are neglected. Because 
the disk surface density fluctuates vigorously with time 
and the migration is stochastic, the accretion rate onto 
the clump also fluctuates. Comparing Eq. 14 with our 
simulations, we found f c <~ 1.7, and 



Our simulations confirm this timescalc estimate early in 
the simulations shortly after clump formation, as shown 
in the right panels of Fig. 8 (the smooth solid curve is 
from Eq. 17). This is consistent with recent papers on 
the same subject by Baruteau et al. (2011) and Michael 
et al. (2011). The upper solid lines in these panels are 
the migration trajectories followed by the clumps in the 
simulations, and the spikes are a stochastic component 
due to the interaction between the clumps and the spiral 
arms (and sometimes other clumps). Our simulations, 
however, suggest that as the clump mass increases due 
to accretion, the migration deviates from the typical type 
I migration time scale given by Eq. 17. This appears to 
occur because the mass of the clump begins to approach 
the local disk mass, such that the inertia of the clump 
plays a role in slowing the migration rate. The type I 
migration timescale given by Eq.17 is obtained using a 
linear theory that assumes that the underlying disc struc- 
ture is unperturbed by the clump. Although the clumps 
are not able to clear gaps in the disk because of the large 
effective viscosity, the disc response to the clump gravity 
is non linear because the clump Hill radii exceed the lo- 
cal disc scale height. We find that the clump migration 
rates can be reasonably-well fit by the following expres- 
sion that accounts for the slowing of migration because 
of the clump inertia 



(15) 



7~mig,fit — ^mig 



l 



in basic agreement with the discussion above. With 
our parameter it is M c = 1.7 x 10~ 5 (R/ 100 AU)- 1 - 25 
(M c /10Mj) 2 ^M Q /yr. Equation 13 shares similar- 
ity with the planet accretion rate from D'Angelo & 
Lubow(2008), but in our cases the Hill radius is larger 
than the disk scale height and the cross section is ~ Th 
instead of ~ rfj . Accretion rates obtained from the simu- 
lations, and estimated using Eq. 15, are shown in the left 
panel of Fig. 8. In our simulations, the accretion rates 
are higher than what Boley et al. 2010 assumed. We 
think this is due to the effect equivalent to gravitational 
focusing. When the flow passes around the planet, it 
is gravitationally focused to the planet so that its cross 
section increases (Fig. 5). Clearly this result is only 
valid if the core radius is smaller than the Hill radius; 
otherwise the clump is tidally limited and accretion is 
suppressed. Further reduction of the Hill radius would 
cause the clump to be tidally disrupted. 

4.3. Clump migration 

In most of our fragmentation cases, the clumps migrate 
inward roughly on the type I timescale initially, given by M Ciini — 7r(A/2) 2 E(i?) 



( M < V 

\2R 2 Y>(R) ) 



(18) 



Fits obtained using the above expression are shown by 
dashed lines in the right panels of Fig. 8. 

The clump represented by the black curve shown in 
the upper panel of Fig. 8 (run R100_3e-5 ) slows down 
its migration when it reaches ^50 AU since it is massive 
enough (~0.3 M©) to open a gap in the disk (see the 
discussion in §5.1). 

4.4. Clump mass 

The clumps accrete mass from the disk during their 
migration as discussed in §4.2. There are several char- 
acteristic masses which mark transitions from differing 
types of evolution, as discussed further in §5. 

The mass limits of clumps in disks undergoing GI can 
be estimated analytically and compared with our simula- 
tion results. The initial clump mass can be estimated as 
the mass contained in the disk patch whose lengthscale 
is equal to that of the most unstable mode (wavelength 
A=2 C 2/GS with Q=l). Thus 



nH 4 M 2 



G 2 S(i?) i? 4 i? 2 E ■ 



(19) 
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In our disk setup, M c , im =4.7x lO^M* (i?/AU)°- 75 (the 
lower dotted curve in the upper left panel of Figs. 9, 
10, & 11). When the clump first forms from the break- 
up of a spiral arm, its shape is elongated due to tidal 
forces from both the star and the rest of spiral arm. This 
elongation together with the shock structure along the 
spiral arms makes the real initial clump mass different 
from that estimated above. 

When the clumps are identifiable and become quasi- 
spherical, separating from the spiral arms, we found they 
are close to that would be implied by Equation 19. But 
some clumps are a factor of 2 less massive (e.g. the upper 
left panel of Fig. 10) 6 . This is slightly higher than more 
detailed clump mass calculations and simulations by Bo- 
ley et al. (2010), suggesting that the initial clump mass 
can be a factor of several smaller than that estimated 
using Equation 19. And our initial clumps are >10 Mj, 
larger than that found by Boley et al. (2010) and Hayficld 
et al. (2011). We attribute this to the strong irradiation 
in our disks requiring larger disk masses before fragmen- 
tation can arise. 

After a clump collapses, its size is smaller than its Hill 
radius, enabling it to accrete from the disk during its or- 
bit around the central star (as discussed in §4.2). If the 
clump accretes all the mass along its orbit, the maximum 
clump mass can be calculated by considering its "isola- 
tion mass" , defined to be the mass contained in an annu- 
lus centered on the orbital radius with half-width equal 
to the clump Hill radius (e.g. Lissauer 1987; Rafikov 
2001): 

and 

(21) 

With our disk parameters, M c ^ so = 
0.023M st (i?/AU) 0375 . This isolation mass is plot- 
ted as the upper dotted line in the upper left panel of 
Figs. 9, 10, and 11. It should be noted that this mass 
is derived assuming that the clump does not migrate 
during accretion. In reality, the clump migrates and is 
therefore able to grow to a mass that is larger than the 
isolation mass. 

In our simulations, the surviving clumps start at a 
mass smaller than the mass estimated by Equation 19 but 
quickly accrete to the isolation mass. This high accre- 
tion rate, and the resulting formation of massive clumps, 
provides a different picture of clump evolution than the 
analytic calculations of Nayakshin (2010a), where accre- 
tion onto the clump was neglected. 

5. CLUMP FATES 

Due to the stochastic nature of the migration and ac- 
cretion, and the different environments in which they 
form, the clumps have different fates in our simulations. 

6 We want to point out that these initial clumps were identified 
by eye, and when we are certain they are gravitationally bound; 
the clumps may be less massive when they initially become bound 
objects. 



In this paper, we have only traced the clumps in the 
marginal fragmentation cases, where 1 or 2 clumps co- 
exist in the disk. When multiple clumps form, clump- 
clump interaction becomes significant and examination 
of this is left for a future study. A total of 13 clumps 
are generated in the marginally fragmenting cases (sum- 
marized in the last column of Table 1 and Table 2) . Of 
these, 3 clumps open gaps and stop migrating (as shown 
in Fig. 9), 4 clumps are tidally destroyed during their 
migration (three are shown in Fig. 10), and 6 clumps 
migrate across the inner boundary (three are shown in 
Fig. 11). 

5.1. Gap opening 

Figure 9 shows the evolution of the clumps that open 
gaps and almost stop migration in the disks at the end 
of our simulations. Results are taken from the simula- 
tions R65_le-4, R100_3e-5 and R200_le-5 listed in Ta- 
ble 1. In the R100_3e-5 case (black curves), we can see 
the trend that the clump's inward migration slows down 
when it opens a gap at ~50 AU. The change in surface 
density during gap formation is shown in Fig. 12. Two 
conditions need to be satisfied for gap opening in a non 
self-gravitating disk. The first condition is that the Hill 
radius of the clump needs to be larger than the disk scale 
height. Using Eq. 19, it can be shown that this condition 
is met when the clump mass M c is smaller than the ini- 
tial clump mass calculated from Eq. 19 by a factor of 12 
(the lower dashed line in the upper left panel of Figs. 9, 
10, and 11). Thus, this condition is satisfied almost all 
the time considering that the initial clump is only a few 
times less massive than the mass given by Eq. 19 . The 
second condition is that M c /M# > 40/i? e (Lin & Pa- 
paloizou 1979) where the effective Reynolds number at 
the clump location is R c = Rv^/ aHc s . Here, is the 
Kcplcrian velocity. For gap opening to occur, the clump 
mass must satisfy the relation 

aH 2 

M c , gap >40M*-^. (22) 

With the irradiation temperature we adopted, H/R ~ 
0.029(i?/AU) - 25 , and M c<gap > 0.034 a M*(i?/ AU) 05 . 
Since a can be as large as 1 (see §6.1), the gap opening 
value of M c is shown as the upper dashed line in the 
upper left panels of Figs. 9, 10, and 11 7 . Fig. 9 shows 
clearly that when the clump mass is above this line, the 
migration slows down significantly. The gap opening is 
illustrated further in Fig. 12, and the gap forms when the 
clump mass M c ~ 0.17M Q , corresponding to the point 
where the clump is close to the critical mass in Eq. 22 and 
stops migrating. Numerous previous studies have showed 
that even in the presence of a gap, gas can accrete onto 
the gap forming object by diffusion of gas through the 
gap (e.g. Bryden et al. 1999; Kley 1999; Nelson et al. 
2000) . It is for this reason that some of the clump masses 
are able to grow beyond the gap forming mass in Fig. 9. 

For the black curve in the top right panel of Fig. 9, the 
clump appears to migrate inward sharply at the end of 
the simulation. By investigating this inward migration 
in detail, we found another clump forms at the edge of 

7 The real gap opening mass can be slightly smaller than that 
shown in these figures since we assume the maximum a=l in these 
plots. 
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the gap, and the interaction between these two chimps 
leads to this inward migration 8 

5.2. Tidal destruction 

Out of 13 clumps in our simulations, 4 are tidally de- 
stroyed during their migration, as illustrated by the up- 
per left panel of Fig. 10). These cases support the clump 
tidal destruction picture suggested by Boley et al. (2010), 
Nayakshin (2010abc), and Vorobyov(2011). As pointed 
out by these papers, since the condition within the clump 
can be very different from the disk, different chemistry 
and solid processing are expected. If the core can be 
tidally destroyed, the products can be released back into 
the protoplanetary disk, affecting planet formation. The 
process of tidal destruction is shown in Fig. 13. Start- 
ing from the upper left panel, the clump is initially well 
within its Hill sphere. When the clump migrates inward, 
its Hill radius decreases until it approaches the clump 
radius (upper right panel). At this stage, due to the 
strong tidal force, the clump is elongated along the ra- 
dial direction (lower left panel). The clump is continu- 
ously stretched and loses mass when it migrates further 
inwards, and finally it dissolves into the spiral arm. 

Whether or not the clump is tidally disrupted depends 
on the competition between cooling/shrinking and mi- 
gration. If the migration is fast, and the clump radius 
becomes larger than its Hill radius (r c > r# ), the clump 
will be tidally destroyed. Assuming r c ~ th when the 
clump first forms, this condition becomes 



or 



-r c < -th = 



Tcool ^ 



r H M c r H R 



3M n 



'mig 



R 



3T m jg/T acc + 1 



(23) 
(24) 



where r cool = -r c /f c due to cooling, r acc = M c /M c , and 

r mig = —R/R, if r c <~ th- The left side is the cooling 
timescale, while the right side is a combination of the 
accretion and migration timescales. As shown above, 
although the migration and cooling timescales are im- 
portant, rapid accretion prevents the clump's tidal de- 
struction since increasing clump mass increases the Hill 
radius. The cooling timescale is given by Equation 47 
below as 



^ = 104 (TcM 



"2-5 / M r 



O.OIMp 



yrs (25) 



if the clump/core is radiative with the polytrope we as- 
sumed (see Appendix C), and 



Wc = 122(^y 18 yrs, 



(26) 



if the clump is convective. During the clump accre- 
tion and contraction, this timescale gets significant larger 
with time. On the other hand, with Eq. 14, the accretion 

8 The interaction between these two clumps can be seen in the 
upper right panel of Fig. 9 where there are periodic oscillation for 
the clump's position in the disk as soon as the second clump forms 
in the disk. The oscillation in the cyan curve is also due to the 
interaction with another clump in the disk. 



timescale is 



7~acc — 



2M 
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M 
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V 100AU 



1.25 



yrs. (27) 



And the migration timescale from Eq. 17 is 

1.75 

r,ni„ = <«4 | ] I I yrs. (28) 



" 784 V 0.0 I.Mq 



Since these three timescales are close to each other, it 
is natural that the clumps may have different fates due 
to the stochastic migration and accretion. Equation 24 
above suggests that rapid migration combined with slow 
accretion can lead to clump tidal destruction. As shown 
in Fig. 14, two similar clumps can form at the same time, 
but the one that undergoes rapid migration is tidally de- 
stroyed, while the one that has a longer time to cool and 
shrink survives. Rapid migration also limits a clump's 
ability to accrete, since it is unable to cool and contract 
during the migration, and ends up more or less filling its 
Hill (Roche) sphere. In the extreme cases this leads to 
tidal disruption as described above. 

5.3. Migration to the inner boundary 

Although 6 out of 13 clumps migrate through the in- 
ner boundary without gap opening and tidal destruction, 
this is simply because we adopt an inner boundary that 
is located at ~ 2.5 to 10 AU, so the ultimate fate of 
these clumps is not known. Some illustrative cases are 
shown in Fig. 11, and are described in Table 1. There 
is one gap opening case and one tidal destruction case 
occuring close to 10 AU, as shown in the Figs. 9 and 
10, showing that both gap opening and tidal destruction 
can occur at small radii. With a smaller inner bound- 
ary, we would clearly expect more clumps to open gaps 
or be tidally destroyed as they migrate into the inner 
disc. It is unclear if any clump can survive migration all 
the way to the central star. The naive timescale argu- 
ments given above seem to suggest that the clump may 
be subject to tidal destruction at smaller radii due to the 
fact that the migration timescale decreases significantly 
at smaller radii with bigger clump mass (M™ 1 ), while 
the cooling timescale increases with bigger clump mass 
(M c 0,5 or independent of clump mass). We note that 
when the clump mass exceeds the local disc mass the 
migration slows down, and the migration time increases 
with clump mass as indicated by Eq.18. Also, because 
the clump central temperature is almost always above 
1500 K when it migrates to the inner disk (see the bot- 
tom panels in Figs. 9, 10, and 11), the clump becomes 
dynamically unstable and should collapse to form a sec- 
ond core (subject to overcoming the angular momentum 
barrier provided by its spin). It will shrink considerably, 
leaving the final fate of the clump unclear. Further ac- 
cretion may lead to gap formation and slowed migration, 
increasing the survival probability of the clump. 

From a statistical point of view, however, we might ex- 
pect a fragmenting disk to form a binary or multiple star 
system. If a less massive clump is tidally destroyed or mi- 
grates to the central star, a new clump can still form in 
the outer disk and migrate inwards. If one clump forms 
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and happens to become massive enough to open a gap, 
it will stay there and survive to form a binary compan- 
ion. The common occurence of close binary star systems 
indicates that inward migration of clumps cannot always 
lead to tidal disruption. 

6. FRAGMENTATION AND CLUMP EVOLUTION: 
ANALYTIC MODELS 

In this section, we will review analytic estimates for 
fragmentation radii and develop a simple model for 
clump evolution. By comparing the numerical results 
with the analytic estimates, we can constrain some fun- 
damental parameters in the analytic model (e.g. we will 
constrain the critical a for the disk fragmentation) and 
study the roles played by different physical processes. 

6.1. Analytic estimates for fragmentation 

The parameter space for Gl-induccd disk fragmenta- 
tion has been explored analytically by several authors 
(Rafikov 2005, 2007, 2010; Levin 2007; Kratter et al. 
2008, 2010a; Cossins et al. 2010; Zhu et al. 2010b). The 
derivations are summarized in Appendix D. 

Analytic estimates for the critical radius of fragmen- 
tation as a function of the mass infall rate are shown as 
curves in Fig. 2 with a c =0.06 or 1, where a c is the criti- 
cal fragmentation viscosity parameter. When the equiv- 
alent viscosity parameter of the gravitationally unstable 
disk is larger than a c , the disk will fragment. Rice et 
al. (2005) has shown that for non-irradiated disks and a 
variety of values of specific heat ratios, a c ~ 0.06. Krat- 
ter & Murray-Clay(2011) suggested that irradiation may 
alter the disk fragmentation criterion. Recently, Rice et 
al. (2011) have done 2-D simulations showing that the 
critical a c decreases with irradiation. 

However, this a c was derived by assuming the disk 
cools at a multiple of the orbital timescale throughout 
fragmentation. With realistic radiative cooling this may 
not be the case because of the trapping of radiative en- 
ergy, as discussed later. 

The trend of these curves can be understood under 
two limits: the disk temperature is dominated by viscous 
heating due to the GI, or it is irradiation dominated. In 
the viscous heating dominated limit, the fragmentation 



radius R* oc a° 3 M 



0.07^1/3^/9 



if oc a c - ivi — ivi t - k r (see Eq. D9). Thus, 

Rf has a very weak dependence on M. This is shown 
in Fig. 2 where Rf is almost vertical at ~ 50 AU when 
the disk accretion rate is high. On the other hand, if the 
opacity depletes by a factor of 1000 (kr is 0.1% of the 
nominal value), Rf decreases by a factor of 5 to ~10 AU. 
In the irradiation dominated limit, Rf oc 

(M/a c )- 4 / 3 Afi- 95 (see Eq. D12). Thus, R f de- 
pends sensitively on M as shown in Figure 2, but is 
independent of the opacity. On the other hand, in 
the envelope irradiation dominated limit (assuming 
T irr =const) , the disk will fragment everywhere as long 
as M > 2cl irr a c /G (Q <~ 1.5 is assumed, corresponding 
to the horizontal dotted line in Fig. 2), below which the 
disk will not fragment at any radius. 

Furthermore, as shown in Fig. 2, Rf is very sensitive to 
the value of a c in the irradiation dominated case (Rf oc 

ai^ 3 ) although it is less sensitive in the viscous heating 
dominated case (Rf oc a°- 3 ). 



Fig. 2 shows that the analytic estimate agrees with sim- 
ulations reasonably well. But it also suggests that a c is 
somewhat larger than 0.06, especially at low infall rates. 
We think this discrepancy is caused by the following: 

1) We assume the infall rates equal the disk accretion 
rates in Eq. D3 and Fig. 2. In reality, in order to con- 
serve angular momentum, some of the infall mass will be 
transported outwards, so the net inward mass flux will be 
smaller than the infall rate. We estimate that all points 
in Fig. 2 should shift down by <0.2 dec if the vertical 
axis is the disk accretion rate. 

2) With irradiation, the disk is gradually built up from 
a gravitational stable state to an unstable state by in- 
fall. This gradual build-up toward GI leads to a transient 
marginal state which avoids the strong density fluctua- 
tions that arise from an initially massive unstable disk 
as considered in most previous simulations. The grad- 
ual build-up of the GI in our simulations makes the disk 
less likely to fragment. This mechanism is similar to that 
discussed by Clarke et al. (2007), who showed that, com- 
pared with the cases starting with rapid cooling, gradu- 
ally increasing the cooling rate suppresses fragmentation. 

3) With realistic cooling instead of orbital cooling, 
radiative trapping makes gravitationally bound clumps 
harder to cool during their collapse, suppressing disk 
fragmentation. For example, the cooling time is defined 
as the ratio of the thermal energy of the clump, E, to 
the clump luminosity, L c , which is given by 



L r = 



If we assume k = KoT a , and E c r^ ~ M c , we have 



1~cool 



4irriaTi- a 



(29) 



(30) 



where cy is the specific heat capacity. If we assume T c r c 
is constant (rolb), the above equation is 



7~cool 



M*C V K T} +a 
Airrffia 



(31) 



Since a — 1.5, t coo i oc T (a more quantitative estimate 
is given in Eq. 47). As the clump collapses and heats 
up, the cooling time gets longer. Clumps with realistic 
cooling have longer cooling times than clumps undergo- 
ing constant orbital cooling in previous simulations used 
to determine a c . Wc: need a smaller t coo i (than that 
used in constant orbital prescription) to start with for 
the clump to continue to collapse. This implies, because 
a oc 1 /t coo i, a larger equivalent a c in the analytic frag- 
mentation estimate. 

Finally, notice that even if a c is the same in the irra- 
diation dominated cases as the viscous dominated cases, 
irradiation still suppresses disk fragmentation (as shown 
by the flattening of the critical R in the irradiation dom- 
inated regime in Eq. D3). This is simply due to irradi- 
ation setting a minimum disk temperature, so that the 
disk needs to be hotter to cool. With the higher inter- 
nal energy and longer cooling time, irradiated disks are 
harder to fragment. From the disk accretion point of 
view, even with the same a c , irradiated disks have larger 
temperature so that the kinematic viscosity v (which 
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is oc T) is larger and thus the critical accretion rate is 
higher. 9 . 

6.2. Analytic models of clump evolution 

Although our simulations are two dimensional, we can 
construct simple analytic models to study the clump evo- 
lution by assuming the clump is spherical. In the follow- 
ing, we will calculate clump evolution using models of 
hydrostatic self-gravitating gaseous spheres to represent 
the first core of the clump, and follow the birthline cal- 
culations as in protostar formation but with migration 
and accretion considered. 

The dynamical timcscalc for clump collapse is much 
shorter than its Kelvin-Helmholtz (cooling) timescale 
(Nayakshin 2010a). Thus we can assume the clump is 
in quasi-hydrostatic equilibrium (the first core) before 
its central temperature reaches ^1500 K when dust sub- 
limates and the opacity drops suddenly (§4.1). Whether 
the first core is convectively stable seems to depend on its 
density, and remains a matter of debate (Wuchtcrl 1993; 
Nayakshin 2010a). In the following we will discuss both 
cases. We want to point out that our radiative core has 
a polytropic index, n e (different from the gas adiabatic 
index), equal to 1.5 and the convective core has a poly- 
tropic index equal to 2.5 (see Appendix C). Boley et al. 
(2010) has found a polytrope with index 2.5 fits Helled 
& Bodcnhcimer (2010) data very well, which is not sur- 
prising considering their core is convective. Notice that 
for the gaseous first cores we are discussing here, poly- 
tropic indices running between 1.5 and 2.5 cover a wide 
parameter space. 10 

As derived from Eq. C2 in Appendix C, the luminosity 
is given by 



647T<7 

3k 



G/j, 



2.5 



0.424 15 Mi- 5 



2.5r, 



-1.5 



(32) 



for a radiative core, and by 



L c , c = 4na (20A8 Ka y 4/5 ^ r^M? . (33) 



for a convective core (Eq. C19). 
With our opacity this can be reduced to 



L r = 2x lO^Lp 



0.01M Q 10AU 



1.5 



and (Eq. C20) 



1.715 x 10 _3 L^ 



\ 2 



0.8 



(34) 



(35) 



0MM Q J \10AUj 

9 The a description for GI assumes locality of the GI (see equa- 
tions Dl and D2), which has been studied by Cossins et al. (2009). 
They have shown that the pattern speed of spiral arms in GI disks 
roughly equals the disk local Keplerian speed, so that the anoma- 
lous term of the energy transport in GI disks is negligible. On the 
other hand, Cossins et al. (2009) and Forgan ct al. (2009) pointed 
out that with higher and higher disk masses, the assumption of 
local energy deposition gradually breaks down since lower mode 
spiral arms become strong. Irradiated disks arc normally massive, 
so our argument is not exact and only acts as a guide. 

10 A polytrope with index 1.5 can also represent the structure 
of a convective core for a monatomic ideal gas sphere with dust 
opacity having temperature dependence power a=1.5 (n e =3-a). 



Thus, the effective temperature is 
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and 
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The clump central temperature is then (Eq. C4 and 
C8) 

O.O1M 



^-^(^(m)" 1 ' (38) 
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The clump loses energy through radiation and shrinks. 
Thus an energy equation is required to derive the clump's 
evolution. The approach we will take is quite similar to 
the protostcllar birthline calculations by Hartmann ct al. 
(1997), but instead of considering energy conservation, 
we consider the conservation of the Jacobi-like energy 
integral (r) in a frame corotating with the clump. We 
have 



r = J d 3 x P Qw 2 + u + w + ^ 



(40) 



where 4>t = — (3/2)f2 2 x 2 is the tidal expansion of the 
effective potential about the clump, u is the internal en- 
ergy per unit mass, and w is the potential energy per 
unit mass. Compared with the total energy calculated 

by 



E 



+ u + w) =K + U + W, (41) 



an additional energy due to the tidal potential is present. 
In the picture of the clump accretion, each of these terms 
sensitively depends on flow structure around the clump, 
which is uncertain. First, due to the tidal force, the ac- 
cretion flow from the outside of the Hill sphere onto the 
clump does not conserve angular momentum with respect 
to the clump unless a special flow geometry is met 11 . 
Thus, the structure within the Hill sphere cannot be de- 
termined by rigorous analytical arguments. Second, even 
if the angular momentum is approximately conserved 
close to the clump and far inside the Hill sphere, ad- 
ditional angular momentum transport mechanism (e.g. 
MM) and the boundary layer physics may be impor- 
tant. For example, if angular momentum within the Hill 
sphere is weak, it is likely to form a rotationally sup- 
ported giant clump. On the other hand, if it is like a 
mini- "protostellar system", the clump is compact and 
slowly rotating. 

For the most simple case, first we study the non- 
rotating hydrostatic core with a polytrope structure hav- 
ing the polytropic index n e , which has the following grav- 

11 In a shearing sheet approximation, the flow needs to be 
symmetric/axi-symmetric in the R direction to conserve angular 
momentum. 
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itational and internal energy 
With 7 = 7/5 for diatomic gas, 



E 



0.5 GMj. 



(43) 



Consider Am amount of gas is added to the clump, 
with the total energy 



AE = -w 



GM c Am 



(44) 



where w represents the coefficient by adding all the fac- 
tors from gravitational energy, thermal energy and the 
kinetic energy brought by the accreted mass (Prialnik 
& Livio 1985). For a cold accretion flow which is in 
Keplerian rotation around the clump nj=0.5, while for 
an accretion disk with boundary layer, zu—1. Then the 
star adjusts to a new polytropic configuration with mass 
M c +Am and radius r c + Ar c . During the time interval 
At of mass addition, the conservation of the Jacobi-likc 
energy L yields 



L c At — 
0.5 



0.5 G{M c + AM c f 



5 - n e 
GM? 



(r c + Ar, 
GM c Am 



ZD — A(f>T 



(45) 



— zu 



GM C M C 



- <PT ■ 



This reduces to 

0.5 GM 2 C . 

L c - Y Lr c + 

5 — n K rj, 

(46) 

At the Hill radius , 4>t is close to the gravitational poten- 
tial with respect to the clump. Deep in the Hill sphere, 
4>t is significantly smaller than u. Thus if the clump 
is very small compared with the Hill sphere, 4>t can be 
ignored, as we will assume in the following. 

As shown in Fig. 6, the rotational energy becomes im- 
portant as the clump accretes a significant amount of disk 
mass, thus it is quite uncertain what is the appropriate 
value for zu. 

If the clump accretion is ignored, the cooling timescale 



is 



Tco ° l - = L = 7 C = 10 HlOAu) UoiMe,' 



t CO o,,c = 122 (^j) 1-8 yrs. 



(47) 
(48) 



In order to form a closed relationship so that we can 
evolve the clump's properties with time, we have made 
several simplifications with attendant large uncertainty. 
The first simplification is to treat the structure of the 
clump as a polytrope. Whether the clump is radiative 
or convective changes the polytropic index significantly 
from ^1.5 to 2.5. To explore this uncertainty, we use 
both radiative and convective models. The second un- 
certainty is the accretion physics/boundary layer around 
the clump. We will explore this by varying zu from 0.5 to 



1. The third simplification is assuming the clump is non- 
rotating, which is not necessarily the case as shown in 
Fig. 6. The rotation of the clump can change our anal- 
ysis significantly, not only because the clump changes 
shape due to rotation, but also the change of (f> T is com- 
parable to u if the clump's radius is comparable to the 
Hill radius. We do not attempt to explore this uncer- 
tainty since there is no simple structure solution for a 
fast rotating clump. Thus, this remains a significant un- 
certainty in our analysis, which the reader should keep 
in mind. 

Equations 15, 32 or 33, 38 or 39, 46 form a closed 
model for clump evolution as long as the initial clump 
properties (M c (t = 0), r c (t = 0)) are given. At time t, 
with known M c (t) and r c (t), the luminosity, L(t), can 
be derived from Eq. 32. Then r c (t + At) can be derived 
based on r c from Eq. 46 with L(t) and M c (Eq. 15). 
M c (t + At) can also be derived with M c from Eq. 15. T c 
from Eq. 38 can be used to derive the opacity and the 
equation of state. At this point, all the quantities for the 
clump at time t + At have been updated. Then with the 
migration velocity given by Eq. 18, the clump's position 
can be updated. By comparing the clump's core size with 
its Hill radius (which depends only on the clump's mass 
and position in the disk), we can determine whether the 
clump can continue accreting or will be tidally destroyed. 

To compare with our simulations, we assume the initial 
clump at 150 AU has mass 0.02 M Q and the clump ra- 
dius is half its Hill radius, as suggested by the simulation 
results. Following the above procedure, the clump evo- 
lution is shown in Fig. 15 with n e =1.5 or 2.5 and w=0.5 
or 1. As shown in the figure, the clump evolution is not 
sensitive to the structure of the clump (n e ), but very 
sensitive to accretion onto the clump (zu). This is due 
to the high accretion rate onto the clump, which domi- 
nates the energy budget in the energy equation above. In 
this sense, the interior structure of the clump may affect 
the clump evolution significantly since a convective core 
can transport angular momentum outwards more eas- 
ily making the clump rotate ridgcdly. Our simple theory 
predicts a higher central temperature compared with our 
simulations shown in Figure 11, which could be due to 
our simple radiative cooling treatment in 2D simulations 
or absence of clump rotation in the analytic models. 

As discussed in §5.2, the interplay between clump cool- 
ing, migration, and accretion determines the clump's evo- 
lution and fate. Cooling allows the clump to shrink and 
inhibits tidal destruction, while inward migration reduces 
the Hill radius and promotes tidal destruction. On the 
other hand, accretion increases the Hill radius and in- 
hibits tidal destruction. Since these three timescales are 
comparable, as discussed in §5.2, stochastic migration 
and accretion may lead to different clump fates. To ex- 
amine the effect of varying parameters in our analytic 
model, we calculated cases with 1/5 and 5 times the mi- 
gration speed. As shown clearly in Fig. 16, rapid migra- 
tion leads to tidal destruction around 30 AU. The clump 
properties here are similar to our simulation results in 
Fig. 10. On the other hand, slower migration allows more 
time to cool, shrink and accrete. Eventually, the clump 
becomes massive enough to open a gap during its migra- 
tion. Clump properties from our simple model also agree 
with those from our simulations of gap-opening clumps 
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in Fig. 9. 

Finally, we want to emphasize that this simple model 
is highly idealized, and the behavior of r c depends on the 
polytropic index, the detailed opacity, and the accretion 
physics around the clump. But it still provides a useful 
tool for exploring different evolutionary scenarios, and 
for comparison with hydrodynamic simulations. 

7. DISCUSSION 

7.1. Planets and multiple star systems 

Recent direct imaging observations have revealed com- 
panions around A type stars (e.g. Fomalhaut, Kalas et 
al. 2008; HR 8799, Marois et al. 2008). Although frag- 
mentation from gravitational instability has been pro- 
posed to explain these giant planets, our results pose 
two potential challenges: 

1. Mass challenge: Including irradiation, the disk needs 
to be quite massive to be gravitationally unstable (M</ > 
O.3M ). Although the initial fragments can have masses 
around a few Jupiter masses, they accrete at high rates 
10 _3 -10 _1 Mj yr _1 , so that after only a few thousand 
years clumps may reach their isolation masses, which are 
in the brown dwarf or low stellar mass regime. 

2. Migration challenge: All clumps migrate inward 
from hundreds of AU to tens of AU over timescales of 10 3 
yr, although the migration has a stochastic component 
due to the interaction with the spiral arms. There are 
at least two distinct fates for these migrating clumps. 
The ones reaching gap opening masses (corresponding to 
subsolar masses in irradiated disks) slow their migration; 
less massive clumps migrate in rapidly and are subject 
to tidal destruction. Generally, only the massive clumps 
survive in our simulations, suggesting that Gl-induced 
fragmentation is best suited for forming massive brown 
dwarfs and close binaries rather than gas-giant planets. 

If binary systems are formed by GI fragmentation in 
disks, the spin axis of these two stars should be parallel 
and perpendicular to the disk plane. Recent observations 
on Kepler 16 systems have suggested the spin axis of the 
binary and the planet around them are all closely aligned 
(Winn et al.(2011)), which support the scenario that they 
are formed in the disk. 

The first cores produced by Gl-induced fragmentation 
can be as large as 10 AU when the core is at 100 AU. 
Its luminosity can be as high as 0.1 L if the core mass 
is 0.1 M or 0.002 L if the core is 10 Mj (see Eq. 34). 
With an effective temperature ~ 100 K (Eq. 36), these 
clumps can be easily detected and resolved by ALMA, as 
illustrated by Fig. 18, but they may be rare because they 
only live for a few thousand years before they collapse to 
the second core stage. 

Boley et al. (2010) and Nayakshin (2010c) have pro- 
posed that dust can settle and agglomerate within frag- 
mented clumps, forming a solid core at the center. Sub- 
sequent tidal disruption of the clump during inward mi- 
gration may then lead to the removal of the gas envelope, 
leaving a surviving rocky body in the disk. Our simu- 
lations suggest the tidal destruction of the envelope can 
happen and may not be rare (4 out of 13 cases) if the disk 
fragments. Whether the solid core can form is beyond the 
scope of this paper (see Cha & Nayakshin 2011). 

Our simulations have several simplifications which re- 
quire further consideration. First, our disks are subject 



to mass infall all the time. If the disk fragments while 
there is no mass replenishment by infall, the clumps at 
the outer edge of the disk may not necessarily migrate 
inwards; in other words, small fragments formed near the 
end of infall may have the best chance of survival at large 
radii. 

A second issue is that we have assumed axisymmetric, 
smooth infall, which may not be the most typical case. 
Strong non-axisymmetric perturbations in the infalling 
material can aid fragmentation (Burkert & Bodcnhcimcr 
1993). We do not explore such large perturbations in 
order to develop well-conditioned numerical simulations. 

We emphasize that the two-dimensional nature of our 
simulations, while enabling us to study clump forma- 
tion and evolution over long timescales, necessarily omits 
three-dimensional effects which could be important, such 
as wave propagation in a stratified disk (Ogilvie & Lubow 
1999). However, the fragmentation occurs at large radii, 
where we expect the disk to be more nearly vertically 
isothermal due to the dominance of irradiation heating. 
Of more concern is the crude treatment of the radia- 
tive cooling for clump evolution, which may differ sig- 
nificantly in three dimensions. Finally, we have limited 
the calculations to a constant value of 7, and have not 
considered possible instabilities related to, for example, 
molecular hydrogen rotational thermalization or dissoci- 
ation. 

Considering our limitations, a comparison with exist- 
ing 3-D simulations is necessary. Boley 2009 has done 
four simulations with similar mass loading as our sim- 
ulations, two of which (SIMB and SIMD) can be com- 
pared with our simulations. SIMB fragments with infall 
rate M ~ 1O~ 4 M yr" 1 at 100 AU, which is consistent 
with our results. And the clump in SIMD accretes from 
6 Mj to 11 Mj supporting the clump accretion scenario. 
Simulations by Kratter et al. (2010a) have shown that 
the disk fragments when a >1 and the clumps are quite 
massive. However they used a piecewise polytropic equa- 
tion of state. Harsono et al.(2011) suggested that infall 
has weakened the disk fragmentation, in other words the 
critical a for fragmentation is large. This seems to be 
consistent with our results. But caution has to be made 
that we put infall manually by adding disk mass and no 
shear between the infall and the disk has been set up. We 
also used radiative cooling instead of orbital cooling as 
Harsono et al. (2011). Hayfield et al.(2011) have found in 
some cases, the clumps accrete at rates ~ lO _5 M yr _1 
and can grow from several Jupiter masses to 40 Mj. This 
is similar to what we have found in our simulations. But 
again, our simple equation of state introduces uncertain- 
ties in the disk dynamics and the thermal structure of 
the clumps. 

Although theoretically the clump evolution and fate 
is complicated by the disk infall history, and clump ac- 
cretion/migration, observations may provide us clues on 
the formation mechanism for giant planets found at large 
distances (e.g. Fomalhaut, HR 8799) by measuring the 
metallicity of these planets (Helled & Bodenheimer 2010, 
2011). 

7.2. Dependence on stellar properties 

Figure 17 shows the analytic estimates for the frag- 
mentation radii for disks around brown dwarfs and A 
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type stars. We chose M*=0.08 M©, L*=0.01 L for the 
brown dwarf, and M*=3 M Q , L*=100 L Q for the A type 
star. The critical fragmentation a c is 1 and minimum 10 
K disk temperature has been applied. 

For viscous heating dominated disks, the fragmenta- 
tion radius varies from 20 AU (brown dwarf) to 80 AU 

1 /3 

(A type star) simply due to Rf oc MJ in Eq. D9. How- 
ever, for irradiation dominated disks, the fragmentation 
radius varies from 30 AU (brown dwarf) to 1000 AU (A 
type star) at 10 _5 M© yr _1 . This sensitive radius depen- 
dence for irradiation dominated disks is due to Rf oc L 1/>2 
in Eq. Dll, and the four orders of magnitude variation in 
luminosity between the A type star and the brown dwarf. 
The lowest fragmentation accretion rate is limited by the 
minimum disk temperature which equals the envelope 
temperature. This analytic estimate is confirmed by the 
R200_3c-5L100 case where the disk did not fragment. 

We want to emphasize that currently found far away 
planets (>50 AU) by direct imaging technique are around 
A type stars (e.g. HR 8799). Due to the stronger irradi- 
ation, the disk needs to be more massive to be gravita- 
tionally unstable and the clumps produced by GI can be 
more massive than our simulation produced. This poses 
a bigger challenge for GI to explain these giant planets. 

These results apply to the formation of fragments for 
infall in a smooth, non-structured pattern, which proba- 
bly applies best to the formation of objects of much lower 
mass than that of the central star. Binary fragmentation 
can be enhanced by non-axisymmetric infall (Burkert & 
Bodenheimer 1993; Kratter et al. 2010a), which are likely 
to dominate the formation of companion stars from tur- 
bulent initial conditions. 

7.3. Disk evolution and outbursts 

Gl-induccd fragmentation has been proposed as a 
mechanism to explain outbursting protostcllar accre- 
tion phenomena, as in FU Orionis and Exor objects 
(Vorobyov & Basu 2005). However, early calculations 
(Vorobyov & Basu 2005, 2006) used a polytropic equa- 
tion of state which can artificially enhance fragmentation 
(Boss 1997; Pickett et al. 2000; Gammie 2001). Recent 
simulations (Voroboyov & Basu 2010) with an improved 
energy equation suggest that Gl-induced fragmentation 
only occurs in a fast rotating massive disk with a low 
background temperature; this is consistent with our re- 
sults that only infall to large disk radii at high infall rates 
results in fragmentation. Thus, whether outbursts can 
be explained by protostellar disks fragmentation sensi- 
tively depends on the initial conditions of prestellar cores, 
which need further observational study. 

On the other hand, while Vorobyov & Basu (2005, 
2006) argued that accretion of fragments explain FU Ori- 
onis outbursts of accretion from disks, we showed (Zhu 
et al. 2008, 2009; see also Armitage et al. 2001, Rice et 
al. 2010, Martin & Lubow 2011) that the GI could trig- 
ger bursts of accretion without clumping by thermally 
triggering the magnetorotational instability (MRI) in the 
innermost disk. We further analyzed FU Ori's SED to 
show that the inner ~ 0.5 AU region of FU Ori during 
outburst is a steady accretion disk and it fits the predic- 
tions of this picture of MRI triggering (Zhu et al. 2007, 
2009). And recently, this ~0.5 AU disk has been con- 
firmed by Keck interferometric observations (Eisner & 



Hillenbrand, 2011). This means that transport in the 
innermost disk is almost certainly dominated by mag- 
netic turbulence rather than gravitational instability, and 
the symmetry of the rotationally-broadened line profiles 
(e.g., Zhu et al. 2009) yields no evidence for the strong 
departure from axisymmetry that would be expected for 
pure clump accretion. 

It is worth emphasizing that we find that some clumps 
can shear out (be tidally destroyed) as they migrate to 
the innermost disk regions if they don't open gaps and 
slow down their migration. Thus, simply finding frag- 
mentation and inward migration is no guarantee that a 
clump will actually accrete onto the central star as an en- 
tity; it is necessary to follow the evolution of the mass in- 
side of 1 AU, which no simulations have done successfully 
yet, and consider the potential if not probable effects of 
magnetic turbulence in addition to GI. It may be, how- 
ever, that the tidal destruction of the clumps produced 
by Gl-induced fragmentation (Boley et al. 2010) or the 
episodic nature of the GI accretion (Boley & Durisen 
2008) can serve as a mechanism to transport mass to AU 
scales (compared with the steady accretion by GI as sug- 
gested by Armitage et al. 2001 and Zhu et al. 2009) to 
trigger the MRI. 

8. CONCLUSIONS 

We have presented two-dimensional hydrodynamic 
simulations of self-gravitating protostellar disks subject 
to axisymmetric, steady mass loading from an inf ailing 
envelope and irradiation from the central star to explore 
the growth of gravitational instability (GI), fragmenta- 
tion, and the fragmented clump evolution. Our results 
can be summarized as follows: 

• We confirm previous analytic estimates that at high 
infall rates, there is a critical disk radius for infall 
<~ 50 AU, beyond which disk fragmentation occurs. 
This critical radius decreases with increasing mass 
accretion rate. 

• At low infall rates, irradiation suppresses fragmen- 
tation and pushes the critical radius to larger val- 
ues. For a solar-mass protostar, this critical radius 
is 200 AU if the infall rate < 1O~ 5 M yr" 1 . 

• With the gradual addition of the mass to the disk 
from the infall, our simulations converge with the 
increasing numerical resolution for both fragmen- 
tation and non-fragmentation cases. 

• The disk critical fragmentation a parameter seems 
to be larger than that found by previous work; we 
suggest this is due to differences in the infall model, 
the thermal history of the disk, and the radiative 
energy trapping of the clumps. 

• For fragmenting disks, the clumps form at a mass 
which is close to (in some cases slightly smaller) 
than the mass in the initial disk contained within 
the typical GI length scale. Rotational support 
becomes increasingly important as the clump col- 
lapses and accretes high angular momentum mate- 
rial from the disk. 

• The clumps accrete from the disk at a rate M c <~ 
4£r^f2, which can be understood by assuming 



Protostcllar disk accretion 



15 



the accretion cross section corresponds to 1-2 Hill 
radii. In our irradiated disks, this is 10~ 5 to 
1CT 4 M yr _1 , suggesting the clump can accrete to 
sub-solar mass on the timescale of thousands of 
years. 

• The clumps migrate inward on the typical type I 
time scale (~ 2 x 10 3 yr in our models) initially, but 
with a stochastic component superposed because 
of interaction with the background spiral waves. 
However, the clump slows down when its mass is 
comparable to the local disk mass, due to its iner- 
tia. A simple fit is given that agrees well with the 
simulation results. 

• When a clump migrates inward its Hill radius de- 
creases. Accretion by the clump, however, coun- 
teracts the effect of the migration and slows the 
shrinking of the Hill radius. 

• There are at least two distinct fates of the clumps, 
depending on the migration speed. If the clump 
migrates slowly with enough time to cool, shrink 
and accrete, it can become massive enough to open 
a gap (M c > O.IMq in our models). This slows 
down its migration (transition to type II migra- 
tion) and may lead to the formation of a binary 
companion (3 out of 13 clumps in our simulation). 
On the other hand, if the clump migrates inward 
quickly, and is not massive enough to open a gap, it 
is subject to tidal destruction (4 out of 13 clumps). 
In the long run, however, the massive, gap opening 
clump should be the final fate for these disks, since 
new clumps can continue forming and be quickly 
tidally destroyed until eventually a massive clump 
forms, opens a gap and remains in the disc for a 
long period of time. 

• A simple analytic clump evolution model with 
clump accretion, migration, and cooling included is 
presented and it confirms that the different clump 
fates depend on the migration and accretion rates. 



• Both rapid clump migration and rapid gas accre- 
tion (leading to clumps much more massive than 
giant planets) pose challenges to the scenario that 
giant planets are formed in situ by GI. And it is 
even more challenging to explain currently found 
giant planets beyond 50 AU around A type stars, 
since, under stronger irradiation, the fragmented 
clumps will be more massive than those produced 
in our simulations. On the other hand, GI frag- 
mentation can provide a formation mechanism for 
close binary systems. 

• The first cores produced by Gl-induccd fragmen- 
tation can be as large as 10 AU when the core is 
at 100 AU. Its luminosity can be as high as 0.1 
L if the core mass is 0.1 M©, or 0.002 L if the 
core mass is 10 Mj. With an effective temperature 
around 100 K, these clumps can be easily detected 
and resolved by ALMA, although these first cores 
may only exist for thousands of years before they 
collapse to form the second core. 

• Although accretion of discrete clumps does not 
agree with observational constraints on FU Ori ob- 
jects, tidal destruction can serve as a mechanism to 
add a stochastic component to mass buildup in the 
inner disk, potentially varying the resulting MRI- 
triggered outburst. 
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APPENDIX 
A: OPACITY TABLE 

The low temperature opacity is important for disk fragmentation. Thus, we have updated our low density opacity 
in Zhu ct al. (2009), including water ice, carbonate and silicate grains. The opacity fits are given in Table 2. 

The dust grain size distribution is chosen to be a power-law with slope 3.5 from 0.005 to 1 /xm. Olivine silicate to 
gas mass ratio is 0.0017, Pyroxene silicate to g itio is 0.0017, graphite to gas mass ratio is 0.0041, and water 

ice to gas mass ratio is 0.0056. We use the dust destruction boundary as in D'Alessio et al. (1992). 

B: BOUNDARY CONDITIONS 

We improve the boundary conditions of the publicly available version of FARGOADSG so that the disk is in pressure 
equilibrium at the boundary. Take the inner boundary as an example. Assuming the disk surface density distribution 
has the same slope with radii in the ghost zone and the computational zone, we extrapolate the density (£) and 
internal energy (E) of the first two computational zones to the ghost zone as 



l °9(Vg,j/Vi,j) = 



logiRi./Ri) 



(logiRg/Ri.)) 



(Bl) 



where r/ij is the density or the internal energy at the grid with radial number i and azimuthal number j and ffi and 
1J2 are the azimuthal averaged value at the first and second radii. 
Then the azimuthal velocities of the ghost zone are calculated by 



iQ = 



E(j - 1) dlogE i?<9$ 
E dlogR + ~dR 



(B2) 
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where Rd<fr/dR is just the GM/i? 9 and dlogE / dlogR is given above. The radial velocity of the ghost zone equals to 
that of the first computation zone only if the radial velocity is inwards, otherwise it is 0. 

With this treatment, we can reproduce the accretion rate of an axisymmetric disk quite well without over-predicting 
the mass accretion rate as the normal free-flow boundary condition leads to. However, the extrapolation may lead to 
a negative density and energy value in the ghost zone if the disk's density and energy fluctuate significantly close to 
the inner boundary. This becomes more severe if a clump forms in the disk and migrates close to the inner boundary. 
Thus for the fragmenting disks, we have also carried simulations with normal outflow boundary condition (Stone & 
Norman 1992). 

C: RADIATIVE CORE LUMINOSITY-MASS RELATION 

Assuming a radiative core (as our 2-D simulations implicitly assume in equation 4), with the uniform energy source 
model and "radiative zero" solution, the clump structure can be approximated by a polytrope P=K' T" e+1 (Chap. 
23.3 of Cox 1968), where n e =— a + 3=1.5 and a is the slope of the opacity on temperature (k=k T q ). Uniform energy 
source model should afford a fairly good approximation of gravitationally contracting stars (Chap. 17 of Cox 1968). 
The structure of the polytropic sphere can be solved by the Lane-Emden equation. A radiative sphere with zero 
temperature boundary condition has luminosity 

L = 64ttct G +4 1 4tt i(-o+3)/n e 

c 3k V n e + l L (n e + l)"^«e+i(_^J^-iJ 

„-a+4 M -a+3 

(Cl) 

where a is the Stefan-Boltzmann constant. For dust opacity a=1.5 in Appendix A. The above equation is 

64^6^ 0.424^ 

With our opacity this can be reduced to 

L c = 2 x 1(T 5 L (#--£-1 15 = 2 x 1(T 3 L ( ^ . (C3) 
\Mq Rq J ° V0.01M o 10AU J K ' 

And the clump central temperature, pressure and density are (notation as in Cox 1968 chap 23.1, and n=n e =1.5, so 
6=3.65375 and -^^=2.71406) 

T C = , * ^ = 1.234 x lO^^^K = 138 *°'™™° K , (C4) 

(n e + l)(-^)i R r c r c /R Q rJlQAU v ' 



with /U=2.4, 



GM c* neR ^ inl7 (M c /M ) 



2 



0.08658 x 10 \ ( \a dyne/cm 2 , (C5) 



Mn e + l)(9') 2 i rt ■ (rc/R Q Y 



Pc = (l/3)(-WiP = 5.99071^3, (C6) 



From equation C4 we can derive 



GM C 1 



»r c T c ^0.537 



(C7) 



With /u=2.4, we can derive GM c /3?r c T c =0.7759. 

On the other hand, for a completely convective core, dlnT / dlnP = (T2 — 1)/T2. With the ideal equation of state for 
diatomic gas (T is high enough that both the rotational and vibrational modes are active), T2 = 7 = 7/5. Assuming 
T2 is constant throughout the core, integrating dlnT/dlnP from the core surface to the center leads to P= K' T™ e+1 
with n e =l/(7 — 1)=2.5. The clump structure is represented by a polytrope P=K' T 3 - 5 which is steeper than our 
radiative case (in stellar cases with monatomic gas, n=1.5, which is the same as our radiative cores). With n e =2.5, 
we have £i=5. 35528 and —£^=2.18720. Thus the central temperature, pressure and density are 

T c = 1.603 x 10 r^/M Q = M./O.OIM^ 

Tc/Rq r c /10AU V 1 



with /U=2.4, 



P c = 0.4935 x 10 iy (Mc ^f rf yne/cm 2 , (C9) 



(C10) 
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which is 



And the K' in the polytrope can be derived using Equation C4 and C5 

K > = (^T g^w 1 i (cn) 

K> ^-5y- 5 3 . 5 5j554^a07626^ 1 
13.65ft 3 - 5 

M 3.5 M 1.5 r 0.5 G 2.5 • ^ °^ 

However the luminosity is determined by the boundary where the convective core becomes radiative. For a fully 
convective core, we assume this occurs at the photosphere where the optical depth t=2/3. The pressure at the 
photosphere is 

/>.,- l'^* 1 ^., ^ -" . (C14) 

Thus 

p eff = _-L^ (C15) 

Then with 

Tpm - nT 7/2 - ' ( C16 ^ 

we can derive 

s /r „\ 3-5/(a+3.5) 
\-l/(o+3.5) / \ -1.5/(o+3.5)^2.5/(a+3.5) 



With our opacity this is 



T eff = (20.48/c )~ n ^ ' r -^/(a+^ )M ^ / (a+^ ) (C1?) 

And the luminosity is 

/r,,\ 14 /(a+3.5) 

L = Aira (20.48 Ko )" 4/(a+3 ' 5) ( ^ ) r (2a+i)/(a+3.5) M io/(a+3.B) _ (C18) 

With a=1.5, 

/r , x 14/5 

L = Ana (20.48« )" 4/5 ( ^ J r c 4 / 5 M c 2 . (C19) 

^ = 1 - 715xlo " 3 ^(o^) 2 (i(S f /) ' 8 - ™ 

which is close to the estimate by Equation C3 at 0.01 M and 10 AU size but with the different dependence on M c 
and r c . 

D: ANALYTIC ESTIMATE FOR THE CRITICAL FRAGMENTATION RADIUS 

The parameter space for GI disk fragmentation has been explored analytically by several authors (Rafikov 2005, 
2007, 2010; Levin 2007; Kratter et al. 2008, 2010a; Cossins et al. 2010; Zhu et al. 2010b). The effects of the irradiation 
and opacity have been discussed in detail in Rafikov (2009) and Cossins et al. (2010). Here we will give a short review 
on these analytic results, so that we can compare them with our numerical simulations. 

To calculate the fragmentation radius at different infall rates analytically we follow Rice (2005) in using the a 
condition for fragmentation instead of using the cooling time condition (Gammie 2001). The a parameter is the 
scaling factor in the Shakura & Sunyaev (1973) treatment of viscosity, v = ac^/fl. 

Rice (2005) has shown that for a variety of values of specific heat ratios (7=5/3,7/5,2), for non-irradiated disks, 
a > a c ~ 0.06 is the condition for fragmentation and it is equivalent to cooling fragmentation condition by Gammie 
(2001). For irradiated disks a c decreases with stronger irradiation (Rice et al. 2011). If the energy balance of the GI 
disk is local, 

aVtclT, ~ 2A C , (Dl) 

where A c represents the cooling of the excess generated energy due to GI heating beyond irradiation heating (defined 
in Eq. 4); then 

t cool = c 2 E/A c ~ -\- • (° 2 ) 
ail 
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In this case if there is a critical fragmentation t coo i, it corresponds to a critical a (Rice et al. 2011). Here, we will argue 
why t coo i = c 2 £/A c is the relevant cooling timescale for fragmentation in irradiated disks, instead of c 2 £/F where F 
is the total outward flux, even though these two definitions degenerate in non-irradiated disks. There are two heating 
sources in the disk: heating due to GI and heating due to the irradiation. The total outward flux F from the disk 
surface also consists of two parts: the part corresponding to the net cooling (A c ) of the energy generated by GI and the 
part balancing the irradiation. Consider a gas parcel in the spiral arm collapses due to its self-gravity, this gas parcel 
gets heated up by the compression. Wether this collapse can continue or not depends on how fast the compressed heat 
can be cooled (as A c in Eq. 4 describes). So that the net cooling ability (A c ) determines the fragmentation instead 
of the total outward flux from the disk. In the extreme case for a locally isothermal disk, A c is infinite, so that t coo i 
defined in equation D2 is 0, and the disk fragments. 

With our definition of t coo i, a and t coo i are directly related and no assumption about the irradiation is needed. In 
the following derivation, we just use the critical a c , as the criterion for the fragmentation in our analytic estimate. 

Since a = MGj2<? s for a Q=1.5 steady accretion disk with accretion rate M and c s , the disk can fragment if 

MG/2c 3 s > a c . (D3) 

Since the disk temperature which determines c s decreases with radii, there is a critical radius (Rf) beyond which 

equation D3 is satisfied and the disk will fragment for a given M. This critical radius has been derived by several 
authors before (Levin 2007; Rafikov 2005, 2010; Clarke 2009; Cossins et al. 2010; Kratter et al. 2010a; Zhu et al. 2010b). 
In order to explore its properties and compare with our numerical simulations, we summarize previous derivations and 
results below briefly. 

Before R/ is determined, the relationship between M and the disk sound speed c s at disk radius R needs to be 
derived first. M, T e ff, and R for a steady accretion disk with zero torque inner boundary are related with 

3GMM ^ 

ei} 8nR 3 a V R ' y ' 



T e ff and T c are related by 

^eff 



rp4 



I^c-^)^, (D5) 



where r = S«;(T c ,p c )/2 and p c = Y./2H. Substituting E with T c by assuming Q = c s VL/-kGY> ~1.5, T c can be solved 
by equating EquationD4 and D5 with the Newton-Raphson method at a given M. Finally, with the derived M-T c , the 
critical fragmentation radii (Rf) can be solved from Equation D3 for any given a c (The results are shown as curves in 
Figure 2). 

The analytical solution of the optically thick limit and the irradiation dominated limit are summarized below: 
With the fragmentation condition assumption MG/2c 3 > a Cl and T c =T C;0 (R/Ro) 7 , the disk will fragment if 

(\ 2/3 7 
s^j ' (D6) 

where c Sj0 is the sound speed at T c , - In the special case where 7=0 (T c =T C fi, e.g. constant irradiation temperature 
at every radius), Eq. D6 is not valid, and the disk will fragment everywhere if and only if M > 2a c c 3 /G 1 independent 
on the disk radius (horizontal part of the curve in Figure 2 at large radii). However, we want to point out that a c 
depends on the irradiation strength (Rice et al. 2011) and in this paper we assume a c is a constant. Thus the key 
parameter to define the fragmentation radii (R/) is the temperature slope 7 which will be derived as following. 

In the viscous heating dominated limit, T 4 =3/16 T 4 ^ Sk, assuming the disk is optically thick which will be justified 
later. Considering the Rossland mean opacity k = n r (T c /T r ) a (T r is some arbitrary temperature to scale the opacity, 
and c Sjr is its corresponding sound speed), T 4 ^=3G]VEM/87rR 3 a, and the Toomre Q parameter of gravitationally 
unstable disk (Q=c s f2/7rGE) is close to some critical value Q c ~1.5, we have 



3 3GM*M c Sir n K r 
16 8?ri?3 crT r 4 ttGQ c 



V{7 ~ 2a) / R \ -9/(7-2o) 



Rd,0 



(D7) 



Please notice that, R^o here could be any value since it is actually canceled out. We keep it here to seperate R from 
the square bracket, but try to maintain the other dimensionless parameters. 

With the dust opacity given in Appendix A, o=1.5, the midplane temperature has slope 7=-9/(7-2a)=-2.25 (T c oc 
i? -2 - 25 ). If we plug Eq. D7 into MGj2(? s > a c , we can derive the disk will fragment if 



/ . x -(14-4o)/27 

R_ I MG \ 

Ro \ 2c 3 s r a c J 



3 3GM*M c s ^K r 
16 8TrR 3 d crT r 4 ttGQc 



2/9 

(D8) 
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Thus, with the opacity given above, we have 




AU as long as the disk accretion rate is high and is viscously heating dominated. The effect of the opacity is also 
apparent: if the opacity depletes by a factor of 1000 (n r is 0.001 of the nominal value), R/ decreases by a factor of 5 
(~10 AU). Finally we will justify why we assume optically thick limit for the viscous heating dominated cases. If we 
plug the opacity law n = K r (T c /T r ) a , and the disk midplane temperature T c =T C! o(R/Rd,o) 7 into Toomre Q parameter 
we can derive the disk is optically thick if R<R 4 , where R t is 



where denotes the corresponding physical quantities at radius R . Assuming T o =220 K at Ro=l AU, and 7=-0.5, 
R t is 80 AU, which is larger than the fragmentation radius of the viscous dominated cases ~50 AU. Thus, it is justified 
that the disk is optically thick at the fragmentation radius. 

For the central star irradiation dominated case, T c — ^T^ rr . Thus with Tf rr — /(i?).L/(47ri?^<7), we have 'y — -1/2 so 
that 




(D10) 




(Dll) 




thus, Rf sensitively depends on M as shown in Figure 2. In this limit, Rf does not depend on the opacity. If T irr =const. 

which is similar to the envelope irradiation case, the disk will fragment everywhere as long as M > 2(? s irr a c /G (where 

c s .i rr is the sound speed at Ti rr , this M is represented as the horizontal dotted line in Figure 2), below which the disk 
will not fragment at any radius. 
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TABLE 1 
Models 



Case name 


Resolution 


Rin 


flout 


Infall R a 


Infall R b 


M in 


Time a 


Fragments ? 


Color in Fig. 8, 9, 10, & 11 






AU 


AU 


AU 


AU 


M Q yr- 1 


10 4 yr 




and comments 


R40_3e-4 


488x512 


2.5 


1000 


25 


40 


3xl0" 4 


0.48 


No 




R40_lc-4 


488x512 


2.5 


1000 


25 


40 


10- 4 


0.4 


No 




R65_3e-4 


432x512 


5 


1000 


50 


65 


3xl0~ 4 


0.1 


Yes 




R65_le-4 


432x512 


5 


1000 


50 


65 


io- 4 


1.24 


Yes/M b 


red, clump A 


R65_3c-5 


432x512 


5 


1000 


50 


65 


3xl0~ 5 


2.04 


No 




R65_lc-5 


432x512 


5 


1000 


50 


65 


10~ 5 


4 


No 




R65_3c-6 


432x512 


5 


1000 


50 


65 


3xl0~ 6 


4 


No 




R100_3e-4 


432x512 


5 


1000 


85 


100 


3xl0~ 4 


0.2 


Yes 




R100_lc-4 


432x512 


5 


1000 


85 


100 


io- 4 


0.67 


Yes 


blue, clump B 


R100_3e-5 


432x512 


5 


1000 


85 


100 


3xl0~ 5 


2.4 


Yes/M 


black, clump C 


R100_lc-5 


432x512 


5 


1000 


85 


100 


10~ 5 


4 


No 




R100_3c-6 


432x512 


5 


1000 


85 


100 


3xl0~ 6 


4 


No 




R200_3c-4 


408x512 


10 


1500 


175 


200 


3xl0~ 4 


0.28 


Yes 




R200_lc-4 


408x512 


10 


1500 


175 


200 


io- 4 


0.47 


Yes 




R200_3e-5 


408x512 


10 


1500 


175 


200 


3X10" 5 


1.2 


Yes 




R200_lc-5 


408x512 


10 


1500 


175 


200 


io- 5 


3.3 


Yes/M 


cyan, clump D,E,F,G 


R200_3e-6 


408x512 


10 


1500 


175 


200 


3xl0~ 6 


20 


No 




Special cases 


R100_le-5I c 


432x512 


5 


1000 


85 


100 


10~ 5 


6.6 


Yes/M 


green, clump H 


R200_3e-5L100 d 


408x512 


10 


1500 


175 


200 


3xl0~ 5 


4 


No 




R200_3e-6noirr ° 


408x512 


10 


1500 


175 


200 


3xl0- 6 


8 


Yes/M 


Magenta, clump I,J,K 


R100_le-5DGlp57 f 


916x1024 


10 


1500 


85 


100 


10~ 5 


4 


No 




R200_le-5Glp57 f 


816x1024 


10 


1500 


175 


200 


10~ 5 




Yes 




R200_lc-5S0p06 « 


816x1024 


10 


1500 


175 


200 


IO" 5 




Yes 




Resolution studies 


R65_le-5D 


816x1024 


10 


1500 


50 


65 


10"° 


3.6 


No 




R65_3e-5D 


816x1024 


10 


1500 


50 


65 


3xl0~ 5 


1.7 


No 




R100_le-5D 


916x1024 


10 


1500 


85 


100 


10~ 5 


3.96 


No 




R100_le-5D2 


704x1024 


20 


1500 


85 


100 


10~ 5 


4 


No 




R100_le-5Q2 


1408x2048 


20 


1500 


85 


100 


10~ 5 


2.65 


No 




R200_le-5 


816x1024 


10 


1500 


175 


200 


10" 5 


3.3 


Yes 





a If fragmentation occurs, this gives time when the first clump migrates to the inner boundary. 
b These are the marginal fragmentation cases where only one or two clumps form in the disk. 
c Decrcasing irradiation after 5xl0 4 yr. 
d Luminosity is 100 Lq from the central star. 
e No irradiation included. 

f The same as the corresponding case in the resolution studies but with 7=1.57 instead of 1.4. 

f The smoothing length is 0.06 disk scale height at 5 AU, compared with other cases with 0.6 disk scale height smoothing length at 5 AU. 





TABLE 2 




Clump Fates 


Gap Opening 


Tidal Destruction Across inner boundary 


a A, C, G 


E, F, I, J A b , B, D, E c , H, K 


3 


4 6 



a The clump label corresponds to the clump label in the last column of Table 1. 

b This clump opens a gap within twice the disk inner radius. The gap opening may be due to the inner boundary condition. Thus we 
count this case as both crossing the inner boundary and gap opening. 

c This clump is tidally destroyed within twice the disk inner radius. This tidal destruction may be due to the inner boundary condition. 
Thus we count this case as both crossing the inner boundary and as being tidally destroyed. 
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TABLE 3 

Fit to Zhu et al. (2007, 2008) opacity but with opacity from water ice, graphite, and silicate dust included 

log 1() T l°gio K comments 

< 0.019 log 10 P + 2.14 1.51og 10 T - 2.48 water ice opacity 

< 0.019 log 10 P + 2.21 -3.53 log 10 T + 0.095 log 10 P + 8.29 water ice evaporation 

< 2.79 1.5 log 10 T — 2.84 metal grain opacity 

< 2.97 -5.83 log 10 T + 17.7 graphite corrosion 

< 0.027 log 10 P + 3.16 2.13 log 10 T - 5.94 grain opacity 

< 0.0281 log 10 P + 3.19 -42.98 log 10 T + 1.312 log 10 P + 135.1 silicate grain evaporation 

< 0.03 log 10 P + 3.28 4.063 log 10 T - 15.013 water vapor 

< 0.00832 log 10 P + 3.41 -18.48 log 10 T + 0.676 log 10 P + 58.93 

< 0.015 log 10 P + 3.7 2.905 log 10 T + 0.498 log 10 P - 13.995 molecular opacities 

< 0.04 log 10 P + 3.91 10.19 log 10 T + 0.382 log 10 P - 40.936 H scattering 

< 0.28 log 10 P + 3.69 -3.36 log 10 T + 0.928 log 10 P + 12.026 bound-free,free-free 

else a —0.48 electron scattering 

a With one additional condition to set the boundary: if log 10 k <3.586 log 10 T -16.85 and log 10 T < 4, log 10 k = 3.586 log 10 T-16.85 
Both C and fortran subroutine can be downloaded at http://http://www.astro. princeton.edu/~zhzhu/opacity.html 
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Fig. 1. — The disk surface density distribution at the end of the simulations (time shown in Table 1) with different infall rates (increasing 
from 3xlO~ 6 M Q yr" 1 on the left to 3xlO- 4 M yr" 1 on the right) and infall radii (65 AU, 100 AU, 200 AU from the top to bottom). The 
black circle labels the Hill radius of the selected clump if the disk fragments. 
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Fig. 2. — The fragmentation radii for gravitationally unstable disks with different infall rates (in unit of Mq/jt). Analytical estimates 
are shown with the critical fragmentation viscosity parameter o c =0.06 (left solid curve) and a c =l (right solid curve). The dotted curve 
represents the fragmentation radii without irradiation (a c =l is assumed). The dashed curves represent the same estimates but with the 
minimum irradiation temperature set to 20 K. Crosses label the disks fragmenting with multiple clumps, while the solid triangle label the 
disk which only has one or two clumps (marginally fragmentation cases). Open squares label the disks which fail to fragment and can 
accrete steadily. 
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Fig. 3. — Clump formation and migration for the R100_le-5I (decreasing irradiation) case compared with the R100_le-5 case (lower right 
panel). For the R100_le-5I case, the upper left panel shows the clump formation and then each consecutive panel shows the clump one 
orbital period later (5.78, 5.89, 6.27, 6.44, 6.53, 6.57, 6.61xl0 4 years after the start of the simulation). After 5 orbits the clump moves 
from 300 AU to 10 AU. The black circle labels the Hill radius of the clump. 




Fig. 4. — The surface density of the clump (run R100_3e-5 ) which just forms (2.2 xlO 4 years after the start of the simulation, left panel) 
and fully developed (2.4xl0 4 years after the start of the simulation, right panel). The black circle labels the Hill radius of the clump and 
the numerical grids have been plotted on top of the density distribution. 
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Fig. 5. — The midplane temperature distribution of the clump as shown in the right panel of Fig. 4 plotted with velocity vectors in the 
frame corotating with the clump center. 
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Fig. 6. — The clump's surface density (dotted curve) and temperature (solid curve) profile along the radial cut across the clump center 
(upper panels) at the two different times corresponding to the left and right panel of Fig. 4. Lower panels show various forces along this 
cut. In the left panels the rotational support is negligible, while in the right panels the rotational support is as important as the thermal 
support. 
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Fig. 7. — Clump's mass (in unit of solar mass), core radius (solid curve), r c /rn (dotted curve), and GM C /5RT C with the clump's position 
in the disk for the R100_3e-5 case. 
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Fig. 8. — Left: The clump mass accretion rate (in unit of MQ/yr) when the clump migrates inwards in the disk. Reading from the top, 
we show cases R100_3e-5 (clump C), R100_le-5I (clump H), R200_le-5 (clump D), R200_3e-6noirr (clump K); Different colors correspond 
to the cases described in Table 1. The dotted lines show the analytic estimates based on Eq. 15. Right: The migration timescale in years 
R/R at each radius for the same cases. The lower solid curves show the type I migration timescales estimated with Eq. 16, while the upper 
dotted curves show the type II migration timescales with a = 1. The dashed curves show our fit using Eq. 18. 
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Fig. 9. — Left: Clump masses and central temperatures as a function of position for cases when a gap is opened and migration stops 
(different colors correspond to different cases in Table 1, and the clump names are given in Table 2). The dotted lines represent the isolation 
mass (Eq. 21) and minimum clump masses (Eq. 19), while the dashed lines represent the gap opening masses (Eq. 22). When the clump 
is close to the inner boundary, its evolution is significantly affected by the inner boundary. Thus we shade the region where the clump 
evolution may not be reliable. The three triangles label when the clump mass is 0.03, 0.08 and 0.17 Mq. The disk surface densities at these 
three times are shown in Fig. 12. Right:The clump radial position in the disk (distance to the star) and the central temperature evolution 
with time for these cases. The black curve is for R100_3e-5 case, the cyan curve is for R200_le-5 case, and the red curve is for R65_le-4 
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Fig. 10. — The same as Fig. 9 but for cases when tidal destruction occurs. These clumps correspond to the clumps labeled in Table 1 
and 2. The cyan curves are for R200_le-5 case, and the magenta curve is for R200_3e-6noirr case. 



Protostellar disk accretion 



31 





Fig. 11. — The same as Fig. 9 but for cases when the clumps migrate to the inner boundary. These clumps correspond to the clumps 
labeled in Table 1 and 2. The magenta curve is for R200_3e-6noirr case, the green curve is for R100_le-5I case, and the blue curve is for 
R10CLle-4 case. 
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Fig. 12. — The disk azimuthal averaged surface density for the R100_3em5 case when the clump mass is 0.03 Mq (solid curve), 0.08 Mq 
(dotted curve), and 0.17 Mq (dashed curve). As clearly shown, the clump stops its migration (corresponding to mass 0.17 Mq in Fig. 9) 
when a gap forms. 




Fig. 14. — Two clumps forming at the same time but with different fates for the R200_3e-6noirr case. From left to right and top to 
bottom, setting the first image at time t = 0, the time slots for subsequent images are 1660 yrs, 3036 yrs, 3084 yrs, 4898 yrs, and 6072 yrs. 
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Fig. 15. — The values of r c /r^, mass (in unit of Mq), and central temperature for a clump as a function of its position in the disk for 
different clump parameters. The purple curves are for the convective cores with n e =2.5, while the cyan curves are for the radiative cores 
with n e = 1.5. The solid curves are calculated with ro=l, while the dotted curves are calculated with ro=0.5. The dashed line in the top 
panel shows when the tidal destruction happens, and the dashed line in the middle plane shows when the gap opening occurs in the disk. 
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Fig. 16. — The values of r c /r^, mass (in unit of Mq), and central temperature for a clump as a function of its position in the disk. The 
clump's evolution is from right to left, while it migrates inwards. The solid black curve is from our simulation with clump K (shown in 
Fig. 8 and 11). The solid purple curve is from our simple analytic model with a convective core assumption and ro=0.5. The dotted red 
curve is from our analytic model but with 5 times the nominal migration speed (Eq. 18), and the clump is tidally destroyed, which can be 
compared with Fig. 10. The dashed blue curve is from our analytic model but with 1/5 of the nominal migration speed, and it open gaps 
in the disk around 30 AU, which can be compared with our simulations in Fig. 9. The dashed line in the top panel shows when the tidal 
destruction happens, and the dashed line in the middle plane shows when the gap opening occurs in the disk. 
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Fig. 17. — Similar to Fig. 2, the fragmentation radi for a typical brown dwarf (0.08 Mq, 0.01 Lq) and an A-type star (3 Mq,100 Lq) 
compared with our standard case (1 Mq, 1 Lq, solid curve). A T = 20 K minimum disk temperature is assumed and a c =l has been used. 
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Fig. 18. — The synthetic ALMA image (with CASA) at 3mm for the R100_3e-5 case, assuming this disk is in Ophiuchus. The disk 
effective temperature at 3mm is approximated with Equation 4 and the central star is assumed to have 5000 K effective temperature. The 
integration is only 1 minute with Full ALMA (resolution 0.1"). Besides the massive clump close to the central star, several other clumps 
at ~ 100 AU scales are also apparent. The tidally induced spiral arms by the clumps have 5 a detection. 



